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We study the linear stability of the flow of a viscous electrically conducting capillary fluid 
on a planar fixed plate in the presence of gravity and a uniform magnetic field, assuming 
that the plate is either a perfect electrical insulator or a perfect conductor. We first con- 
firm that the Squire transformation for magnetohydrodynamics is compatible with the 
stress and insulating boundary conditions at the free surface, but argue that unless the 
flow is driven at fixed Galilei and capillary numbers, respectively parameterising grav- 
ity and surface tension, the critical mode is not necessarily two-dimensional. We then 
investigate numerically how a flow-normal magnetic field, and the associated Hartmann 
steady state, affect the soft and hard instability modes of free-surface flow, working in 
the low-magnetic-Prandtl-number regime of laboratory fluids (Pm ^ 1CP 4 ). Because it is 
a critical-layer instability (moderately modified by the presence of the free surface), the 
hard mode is found to exhibit similar behaviour to the even unstable mode in channel 
Hartmann flow, in terms of both the weak influence of Pm on its neutral-stability curve, 
and the dependence of its critical Reynolds number Re c on the Hartmann number Ha. 
In contrast, the structure of the soft mode's growth-rate contours in the (Re, a) plane, 
where a is the wavenumber, differs markedly between problems with small, but nonzero, 
Pm, and their counterparts in the inductionless limit. As derived from large- wavelength 
approximations, and confirmed numerically, the soft mode's critical Reynolds number 
grows exponentially with Ha in inductionless problems. However, when Pm is nonzero 
the Lorentz force originating from the steady-state current leads to a modification of 
Re c (Ha) to either a sublinearly increasing, or decreasing function of Ha, respectively for 
problems with insulating and conducting walls. In the former, we also observe pairs of 
counter-propagating Alfven waves, the upstream-propagating wave undergoing an insta- 
bility driven by energy transfered from the steady-state shear to both of the velocity and 
magnetic degrees of freedom. Movies are available with the online version of the paper. 



1. Introduction 

Free-surface shear magnetohydrodynamic (MHD) flows arise in a variety of industrial 



and astrophysical c ontexts, including liquid-metal bla nkets in fusion devic es (jAbdou et al 



l200ll lBuhlerll2007T) . liquid- metal fo rced-flow targets (IShannon et al\\ 19981) . plasm a oceans 
on white dwarfs and neutron stars ( Bildsten fc Cutler 1995t lAlexakis et aZj|2004) 1 and ac- 
cretion discs around stellar remnants (jRiidiger et all 1999 ; Balbus fc Henrill2007t l . Flows 
of this type typically take place at high Reynolds number Re > 10 4 and within strong 



2 



D. Giannakis, R. Rosner and P. F. Fischer 



background magnetic fields (Ha > 10 2 , where Ha is the Hartmann number). Broadly 
speaking, one is interested in characterising their stability properties either because of 
engineering requirements (e.g. in a fusion-device blanket), or a free-surface instability 
may be involved in the observed phenomena (such as classical novae and neutron star 
X-ray bursts). 

When the magnetic Prandtl number Pm of the working fluid is small, the effect of 
an exte rnal magnetic field is g enerally known to be stabilising and weakly dependent 
on Pm ( Miiller &: Biihler 2001 , and references therein). In particular. iTakashima ( 19961) 
has numerically studied the stability of plane Poiscuille flow modified by a flow-normal 
magnetic field, hereafter called channel Hartmann flow, under insulating boundary condi- 
tions, and has determined that the critical Reynolds number Re c for instability increases 
monotonically with the Hartmann number Ha. Moreover, for Pm ^ 10 (an interval 
encompassing all known laboratory fluids) Re c was found to experience a mild, O(10 -3 ), 
decrease compared to its value in the inductionless limit Pm \ 0. 

In the present work, we pursue a similar stability analysis for free-surface Hartmann 
flow, i.e. the parallel steady flow of a viscous electrically conducting capillary fluid on an 
inclined planar surface (assumed to be either a perfect electrical conductor, or a perfect 
insulator), subject to gravity and a flow-normal external magnetic field. Our main result 
is that for large wavelengths (a < 1, where a is the modal wavenumber), the spectrum 
of free-surface Hartmann flow contains two types of normal modes, neither of which is 
present in channel problems, and whose stability depends strongly on Pm, even in the 
Pm = 0(1O~ 5 ) regime of liquid metals. 

The first of these modes is related to the unstable gravity wave, often times referred to 
as the soft instability mode, encountered in fre e-surface Poiseuille flow (|Yihlll963l 1969; 



Lam fc Bavazitoglu 19861 : Florvan et al. 1987 ). In inductionless problems, the gravity 



wave becomes stabilised when Ha is increased, but in the strong-field limit its growth rate 
r does not obey the characteristic \T\ oc Ha 2 Lorentz damping of the shear modes in the 
spectrum. Instead, it transitions to an asymptotically neutral phase, where the Lorentz 
force nearly balances gravity, and the decay rate \r\ cx Ha" 2 decreases towards zero. 
In problems with nonzero, yet small, Pm, Lorentz forces originating from the Hartmann 
current profile are found to be sufficient to alter the near-equilibrium attained in the 
inductionless limit, leading to high sensitivity of the mode's stability contours in the 
(Re, a) plane to the magnetic Prandtl number. 

The second type of modes in question is travelling Alfven waves, which we have only 
encountered in flows with an insulating lower wall. When the background fluid is at 
rest, these modes appear at sufficiently large Hartmann numbers as a pair of counter- 
propagating waves, whose phase speed and decay rate increase linearly with Ha. Their 
kinetic and magnetic energies are nearly equal, but in the large- wavelength cases studied 
here the majority of the magnetic energy is carried by the magnetic field penetrating 
into the region exterior to the fluid. When a steady-state flow is established, the up- 
stream propagating Alfven mode becomes unstable due to positive Reynolds and Maxwell 
stresses as the Alfven number Al is increased. 

Aside from instabilities associated with gravity and Alfven waves, free-surface flow 
also exhibits an instability of the critic al-layer type (modified by the presence of th e 
free surface), called the hard instability ( Lirj|l967l ; iDe Bruinlll974t iFlorvan et aZ.lll987 ). 
Sharing a common origin with the even unstab le mode in channel Hartmann flow (jLock 
ll955HPotter" fc Kutche vl ll973t lTakashimafl 996) , the critical parameters of the hard mode 
at small Pm are close to the corresponding ones in the inductionless limit. However, 
in light of the presence of gravity and Alfven waves, our analysis suggests that the 
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Figure 1. Geometry of free-surface (left) and channel (right) Hartmann flow. The steady-state 
velocity and induced magnetic field profiles are U(z) and B(z), respectively (see A2A\ . The y 
axis is directed into the plane of the paper. 



inductionless approximation must be used with caution when dealing with free-surface 
MHD. 

The plan of this paper is as follows. In fj2] we formulate the governing equations and 
boundary conditions of our stability problems, and discuss the validity of the Squire 
transformation. In iJ3]wc derive an energy-conservation law for temporal normal modes in 
free-surface MHD. We present our results in and conclude in $5] Appendix |A~1 contains 
a discussion of large-wavelength (a \ 0) perturbation theory. Movies illustrating the 
behaviour of the modal eigenvalues on the complex plane as Ha or Pm are varied are 
available with the online version of the paper. 



2. Problem formulation 

2.1. Geometrical configuration 

Using x, y, and z to respectively denote the streamwise, spanwise, and flow-normal 
coordinates, oftentimes collected in the position vector r := (x,y, z), and t to denote 
time, we consider the flow geometries shown in figure [TJ In free-surface problems the 
lower planar surface z = — I is inclined at an angle 9 with respect to the horizontal, and 
the upper fluid boundary has oscillation amplitude z — a(x, y, t), with a = in the steady- 
state, while in channel problems the flow takes place between two fixed parallel plates 
located at z = ±Z. In both cases, the working fluid is incompressible, and has density 
p, dynamic viscosity /i, and electrical conductivity A. Additionally, the free surface has 
surface tension a, and is also acted upon by a gravitational field g :— g(sin(6)x—cos(6)z). 
The fixed plates are treated either as perfect electrical insulators or perfect electrical 
conductors. 

For future convenience we introduce the function A(r, t) := z— a(x, y, t), which vanishes 
on the free surface, and leads, through its gradient, to the expression 

n := VA/\\VA\\ = (-d x a, -d v a, 1) + 0(a 2 ) (2.1) 

for the free-surface outward unit normal (for our purposes it suffices to work at linear 
order in a). Moreover, we choose := (1,0, d x a) and t^ y ' := (0,1, d x a) as mutually 
orthogonal unit vectors tangent to the free surface [tS x ^ -n — v- v > -n = •t^ v ' = 0(a 2 )). 
The divergence of n is equal to twice the mean surface curvature k, for which we compute 



2k = V • n = -{d 2 x a + d 2 a) + 0{a 2 ). 



(2.2) 
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2.2. Governing equations 



Our s tarting point is the equations for incompressible resistive MHD (e.g. lMiiller fc Biihler 



d t U + U-VU = -W + F + Re' 1 MA, 
.-U 



200lh . 

F :=RmJxB, J := Rm~ L V X B = £ + U X B, 

d t B + U-VB = B-VU + Rm- 1 AB, 

V-U = 0, V-S = 0, 

obeyed by the velocity field U(r,t), the hydrodynamic pressure V(r,t), the Lorentz 
force J-{r,t), the current J"(r,t), and the magnetic and electric fields in the interior of 
the fluid, respectively B(r,t) and £(r,t). Here velocity has been non-dimensionalised 
by its steady-state value at z = 0, £/*, and the characteristic values for the remaining 
dynamical variables are P* := pU 2 , B* := (poP^^U*, -E* := t/*-B*, and J* := \E*, 
where each *-subscripted symbol denotes the characteristic value for the corresponding 
variable in script type. Choosing I as the characteristic length (for both free-surface 
and channel problems), the resulting hydrodynamic and magnetic Reynolds numbers 
are Re := UJ/v and Rm :— UJ/r], where v ;= p/p and 77 := l/(/i A) are the viscous 
and magnetic diffusivities. In the following, we frequently substitute for Rm using the 
magnetic Prandtl number Pm := v jr\ = Rm/ Re. 
We consider solutions of the form 

U(r,t) = U(z)+u(r,t), V '(r,t) = P< \x,z) + p' \r,t), 

B(r,t) = B(z) + b(r,t), J(r,t) = J(z) + j(r,t), S(r,t) = E{z) + e(r,t), 

consisting of steady-state components and linear perturbations, respectively denoted by 
uppercase and lowercase symbols. The steady-state flow U(z) := (J7(z),0,0) is assumed 
to be streamwise-invariant and unidirectional, and to take place within a uniform, ex- 
ternally applied magnetic field B' :— (B' x , B' y , B' z ), which, for the time being, is allowed 
to be of arbitrary direction. The applied field permeates the fluid, and assuming its 
components perpendicular to U are non-zero, the associated induced current generates 
a secondary internal magnetic field J := (I x (z), I y (z), 0). Thus, in the interior of the 
fluid the steady-state magnetic field is B := B' + I. We remark that since pressure 
affects the dynamics only through its gradient, we have allowed P' to depend on the 
streamwise coordinate x, but, in light of the streamwise-invariance of the steady state, 
that dependence can be at most linear. Moreover, the flow-normal component of the 
induced magnetic field has been set to zero in order to meet the divergence-free con- 
dition V • I = (a constant nonzero I z can be absorbed in B' z ). We also note that 
with our choice of characteristic magnetic field, tl and B are naturally additive. In par- 
ticular, using B' and B^ to respectively denote a unit vector in the direction of B' 
and the magnitude of the dimcnsionful external magnetic field, B can be expressed as 
B' = B ' I A /. where Al := U^n^p) 1 / 2 / B'^ is the Alfven number of the flow (see e.g. 
ISherclifj | 19651 ) for an overview of dimensionless groups in MHD). An alternative option 
for non-dimensionalisation, frequently encountered in the literat ure for Hartmann flow 
(e.g. lTakashimalll996l Il998t iMuller fc Biihlerll200ll: lBuhler|[2007h . is to set B, = B 1 ^ in 



which case the resulting dimensionless magnetic field B is related to ours via B — AlB. 
In the ensuing analysis, we mainly parameterise the background magnetic field strength 
by means of the Hartmann number Ha := B'J(X/ /1) 1 / 2 = {ReRm) 1 / 2 / Al, where Ha 2 
measures the ratio of Lorentz to viscous stresses, rather than Al. 

In problems with insulating boundaries, a further dynamical variable is the magnetic 
field B'(r, t) := B' + b'(r, t) in the region exterior to the fluid. As follows from Ampere's 
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law, b' is expressible as the gradient of the magnetic potential ip(r,t), which, in light of 
the solenoidal condition V • 6 = 0, obeys Laplace's equation; i.e. 

b' = -W, with Aip = 0. (2.5 a, b) 

The equations governing the steady state and the perturbations follow by substitut- 
ing (|2.4|) into (|2.3p . and neglecting quadratic terms in the perturbed fields. Using D to 
denote differentiation with respect to z, the nonzero components of the time-independent 
equations read 

Re~ 1 D 2 U + B' z DI x -d x P = 0, Rm^ 1 D 2 I x + B' z D U = 0, DP = 0, (2.6 o-c) 

B' z D I y = 0, D 2 I y = 0, (2.6 d, e) 

where P(x,z) := P'(x,z) + B(z) • B(z)/2 is the total steady-state pressure, consisting 
of hydro dynamic and magnetic contributions. As for the perturbations, we obtain 

d t u + Ud x u + u z D Ux = -Vp' + f + Re~ 1 Au, (2.7 a) 

f := Rm(j X B + J X 6), j = Rm' 1 V X b = e - Ub z y + u X B, (2.7 b, c) 

<9 t b + J79 x b + u z D B = B • Vu + b z D J7cc + Rm^Ab, (2.7 d) 

V-« = 0, V-b = (2.7 e,/) 

where / is the Lorentz force acting on the perturbed velocity field. 

From the linearised Ampere and Ohm laws (12. 71 c). one can make the heuristic esti- 
mate ||6||/||u|| = O^PmHa) 1 / 2 ) (see HuntJ ([l966) for similar scaling arguments), where 



|| • || denotes suitable norms for the velocity and magnetic field perturbations. This sug- 
gests that as Pm \ with Re and Ha fixed, i.e. in the so-called inductionless limit 



( Miiller fc Buhlcr 200l|), b is negligible and, in consequence, electromagnetic forces only 
arise from currents generated by the perturbed electric field e, and from currents induced 
by the perturbed fluid motions within the steady-state magnetic field. Moreover, as fol- 
lows from Faraday's law, V X e = —dtb w 0, the electric field can be determined from 
the gradient of a potential cf>(r,t). These observations suggest that for sufficiently small 
Pm, (|2.7l 6-d. f) can be replaced by the approximate relations 

/ = Ha 2 Re~ 1 (e + u X B') X B' , e = -V<p, A(j) = 0. (2.8 a-c) 

If, in addition, the flow is two-dimensional, <f> must be set to zero (as <fi ^ would 
give rise to spanwise forces), resulting in a signi ficant reducti on of the analytical and 
computational complexity of the stability problem ( StuartJ[l954|) . The fact that all known 



conducting laboratory fluids have small magnetic Prandtl numbers (Pm < 10 -5 ) has led 
to a widespread adoption of the inductionless approximation. However, the small- 1 1 b 1 1 
assumption is not guaranteed to hold a priori, and the full problem must be solved to 



confirm that the scheme is valid in the parameter regime of interest (jTakashimal 1 1 9961 ) . 

2.3. Boundary conditions 

The governing equations presented in the preceding section must be solved subject to 
appropriate initial and boundary conditions. In the temporal stability analysis that fol- 
lows initial conditions, as well as periodic boundary conditions on the streamwise and 
spanwise domain boundaries, are implicitly assumed. However, care must be taken in the 
choice of bo undary con ditio ns on the non-periodic bo undaries, as this has led to errors 
in t he past (TLinl (Il967l) and lPotter fc Kutchevl (|l973f ). as indicated by |Pe Brum] (|l974f) 



and Takashimal (|l996l ). respectively). 



Let z w collectively denote the flow-normal wall coordinates (in the dimensionless rep- 
resentation, z w := { — 1} for free-surface problems and z w := {±1} for channel problems). 
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In insulating-wall problems we assume that no surface charges and surface currents are 
present at the fluid-wall interfac e. Then, in ac cordance with Maxwell's equations and 
charge-current conservation (e.g. ISherchflll965T ). we set 

B\ z=Zw = B'\ z=Zw , z ■ J\ z=Zvj = 0, (2.9) 

which leads to the boundary conditions 

l x (z w )=0, l v (z w )=0 (2.10 a, 6) 

and 

(& + VV0|*=.„ =0, z-V Xb\ z=Zw = (d x b v -d y b x )\ z=Zw =Q, (2.11 a, b) 

respectively for the steady-state fields and the perturbations. If, on the other hand, the 
wall is conducting, the tangential electric-field components are required to vanish at the 
boundary, and the wall-normal magnetic field z • 6 is set to the externally imposed value 
B' z , giving 

DI x (z w )=0, DI y (z w ) = 0, (2.12 a, b) 

and 

b z \ z = Zw =0, x- V X b\ z=z . w = (d y b z -Db y )\ z = Zw =0, 

y V x b\ z=Zm = (D b x - d x b z )\ z=Zw =0, ( ' } 

where we have used (|2.7l c) to substitute for the electric field e in terms of b. Turning now 
to the free surface, we assume throughout that the exterior region z > a is electrically 
insulating. Then, on the basis of similar electrodynamic arguments as those used to write 
down (|2.9p . we demand 



tW-(B-B%= a =tM-(B-B%= a =n-(B-B')\ z = a =0, n-VxJ\ z=a = Q. (2.14) 

In order to evaluate the expressions above for the perturbed quantities, we analytically 
continue the induced magnetic field I in the region < z ^ a. A Taylor expansion to 
linear order in a then yields the boundary conditions 

4(0) =0, I y (0) = 0, (2.15 a, b) 

and 

(d x ^ + b x )\ z = Q + DB x (0)d x a = 0, (d v ip + b y )\ z=0 + T) B y (0)d y a = 0, 
(D V + b z )\ z=0 = 0, (d x b y - d y b x )\ z=0 + D B y (0)d x a — DB X {0)d y a = 0, U " ° ' 

which now involve the free-surface amplitude (cf. (|2.11|) V The requirement that, when- 
ever present, the external magnetic field perturbations vanish at infinity completes the 
specification of boundary conditions for the magnetic field. 

Regarding the velocity field, we impose, as usual, no-slip conditions 

U(z w ) = 0, u\ z=Zm =0, (2.17 a, b) 

at the solid walls, and consider the kine matics and stres s balance to establish boundary 
conditions at the free surface (see e.g. iBatchelorl I1967L for a discussion on interfacial 
dynamics). First, noting that the free surface is, by definition, a material surface, leads 
to the kinematic boundary condition, dA/dt := (dt + U\ z=a • V)A = 0, which, upon 
linearisation, becomes 

d t a + U(0)d x a = u z \ z = . (2.18) 
In order to formulate the stress conditions, we introduce the stress tensors in the fluid 
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and exterior domains, whose components in the (x, y, z) coordinate system are given by 

%i ■■= -V&V + + Re^&j, % ■= -(* + "PB')Sa + BiB'j, (2.19) 

respectively. Here V := V' + B • 8/2 is the resultant of the hydrodynamic and magnetic 
pressures, 5jj := diUj + djUi are the components of the rate-of-strain tensor, Vb 1 {r , t) := 
B' • S'/2 is the external magnetic pressure, and <P(r) :— (— xsinO + zcos8)/Fr 2 is the 
gravitational potential, expressed in terms of the Froude number Fr :— U r /(gl) 1 ^ 2 . Us- 
ing We := plU 2 /cr to denote the Weber number, the free-surface curvature k (|2.2p . in 
conjunction with surface tension, introduces a discontinuity E := 2k/ We in the normal 
stress, such that 

n j {T( j -T ij )\ z = a = En h (2.20) 

where are the components of the normal vector n (|2.ip . and summation is assumed 
over repeated indices. Forming the contraction of (|2.20p with the orthonormal vectors 
n, v- x \ and then leads to three stress conditions that we enforce at the free surface, 
namely the normal-stress boundary condition 

n i n j {Tl j -T ij )\ z = a =E, (2.21) 

and the shear-stress conditions 

tfn^ - % 3 )\ z=a = t^n.iXj - % 3 )\ z=a = 0. (2.22) 

Evaluating (|2.21l) and (|2.22|) to linear order in the perturbed quantities, and eliminating 
b'i\ z =a using (|2.16l) yields 

P(x,0) = -sm{e)Fr- 2 x + B{0) -B(0)/2, D U(0) = 0, (2.23 a, b) 

and 

(P-PB- 2Re- 1 Du z )\ z=0 = {cos(9)Fr- 2 + JB(0) •DB(0))a 

- 2RC- 1 D U(0)d x a - We^^a + d 2 y a), (2.24) 

(d x u z + Dm x )| 2=0 = -D 2 J7(0)a, (d y u z +Du v )\ z=0 = 0, 

where ps '■= B • b is the internal magnetic-pressure perturbation. 

Besides We and Fr, alternative dimensionless groups for the capillary and normal 
gravitational forces are the capillary number Ca :— fiU*/a = We/Re and the Galilei 
number Ga := gl 3 cos0/v 2 = Re 2 cos(6) / Fr 2 . As we will see in tj2.5l unlike We, Fr 
and 9, the parameters Ga and Ca are invariant under the Squire transformation from 
three-dimensional to two-dimensional normal modes, and for this reason we have opted 
to perform the calculations in S|4] working in the latter representation. 

2.4. Steady-state configuration 

In the present linear-stability analysis we employ the physically motivated Hartmann 
velocity and magnetic field profiles, which are solutions to (|2.6l a-c) (in anticipation of 
the Squire transformation in £|2.5[ we do not explicitly consider the spanwise induced 
magnetic field I y ). For convenience we make the substitutions I x (z) = B' z RmB(z) = 
Pm 1/2 H z B(z) and, taking into account (JM]c), P = -IIx+P , where H z := (ReRm) 1 / 2 B' z , 
-ZT, and Pq are respectively a Hartmann number defined in terms of the flow-normal com- 
ponent of the applied magnetic field, a streamwise pressure gradient, and an unimportant 
constant. Equations (|2.6l a. b) then become 



D 2 U + H 2 D B + ReM = 0, D 2 B + D U = 0, 



(2.25) 
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the general solution of which can be expressed as 

TT(v s r r smh{H z z) cosh(g 2 ) - cosh(g 2 z) 

U(z) = C + Cx ^ + C 2 cQsHHz) _ i , (2.26 a) 

nf \ v _l ^ iP 1 ~ cosh ( H z z ) | ^ sinh(ff z z) - sinh(g z )z 

= + K x z + C X + C 2 Hz{cosh{Hz) _ 1} , (2-26 6) 

where Co, Ci, C2, -Ko, and K\ are constants such that 

n = (C 2 H z coth{H z /2) - H%Ki)/Re. (2.27) 

We remark that the terms proportional to C\ correspond to antisymmetric velocity and 
symmetric magnetic field solutions with respect to z, whereas the constant C 2 gives 
rise to symmetric velocity and antisymmetric magnetic field profiles. Moreover, the term 
proportional to K\ in (|2.26l b) can be interpreted as the magnetic field due to a uni- 
form current of magnitude Ki/Rm, driven in the spanwise direction. For values of the 
Hartmann number approaching zero (|2.26[) . becomes 

U(z) = C + Ciz + C 2 (l - z 2 ) + 0(H 2 Z ), (2.28 a) 

B(z) = K + K lZ - C\z 2 /2 - C 2 z{l - z 2 )/Z + 0(H 2 ), (2.28 6) 

and (|2.27j) reduces to 77 = 2C 2 / Re. As expected, in this non-MHD limit U is a quadratic 
function of z and, even though B(z) is nonzero, the streamwise induced field I x = 
Pm 1 H Z B vanishes. 

In both of the free-surface and channel problems considered here, the choice of veloc- 
ity normalisation (U(0) = 1) and boundary conditions (equations (12.171 a) and (I2.23I M) 
implies that Cq = C\ = 0, and C 2 = 1. Moreover, if the walls are insulating the 
constants Kq and K\ vanish due to (|2. IOI q') and (12.151 a 1 ). In conducting- wall problems 
Kq is again set to zero, either because of (|2.15l a). or, in channel problems, by con- 
vention (a nonzero K can be absorbed in the applied magnetic field B' x ). However, 
Ki = (smh(H z ) — H z cosh(H z ))/(H z (cosh(H z ) — 1)) is in this case non-vanishing, due 
to (HHUa). 

Figure [2] illustrates the functional form of these two classes of velocity and mag- 
netic field profiles for Hartmann numbers in the range 0-20. Compared to the parabolic 
(Poiseuille) profile in non-MHD flows, Hartmann velocity profiles are characterised by a 
flat core region and exponential boundary layers of thickness 0(1/H Z ) near the no-slip 
walls. Moreover, the mean steady-state velocity, given by 

J -1 cosh(i/ z ) - 1 

increases monotonically from its non-MHD (H z = 0) value of 2/3 to unity as H z — » 00. 
In problems with conducting walls, |-B(z)| peaks at \z\ = 1, and its gradient, which is 
proportional to the spanwise induced current J y = i?m _1 D/ 2; = Hz(ReRm)^ 1 / 2 D B, 
attains its maximum magnitude at z = 0, where |D£>| = 1. Also, for large H z the 
current is close to its maximal value throughout the core. On the other hand, if the wall 
is insulating, the current distribution becomes concentrated over the Hartmann layer 
as H z grows, with the magnetic-profile gradient reaching its maximum absolute value 
|D5(±1)| = (cosh( J ff z )-sinh(^)/ J ff 2 )/(cosh(^)-l) = 0(1) at the walls, while at the 
core, where | D B\ = 0(1/ H z ), the current tends to a /^-independent value. 

Before proceeding, it is worthwhile to note a prominent qualitative difference between 
MHD and non-MHD velocity profiles, which concerns the existence of inflection points. 
Because the non-MHD velocity profile (|2.28l a) is quadratic in z, it has constant second 
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Figure 2. Steady-state velocity (left) and magnetic field profile (right) for Hartmann flow with 
H z 6 {0, fO, 20}. In the right-hand panel, the solid (dashed) lines correspond to insulating 
(conducting) boundary conditions at z = — 1. 



derivative. On the other hand, it is possible to show that for suitable choices of the 
constants C\ and C2, and for sufficiently large H z , the MHD velocity profile (12. 261 a) 
possesses an inflection point, so that an inviscid instability can potentially exist. In all 
of the flows studied here, the choice of boundary conditions completely suppresses the 
antisymmetric component of U (C± = 0), and eliminates the possibility of inflectional 
instabilities. However, one can imagine situations (e.g. a sheared free surface) where the 
conditions for the inflection point to exist are satisfied. Whether the non-ideal MHD flow 
develops in practice an instability mode of inviscid origin would be an interesting topic 
to investigate in the future. 

2.5. Three-dimensional normal modes and the associated Squire transformation 
In the three-dimensional tcmporal-normal-mode analysis we work with the Ansatz 

u = Re((u.Az),u y (z),M^)y {ax+0y)+lt ), P = M(p(z> i{ax+M+7 % 
b = Rc((b x (z), b y (z)Mz))e i{ax+0y)+ ^), 4> = ReWzy^+Py^ 1 ), 

where p := p' + B • b is the linear perturbation of the pressure field V, Ui, p, bi, and 
ip are complex functions of z, a and fj are respectively the streamwise and 
spanwise wavenumbers, and 7 £ C is the complex growth rate. In channel problems with 
conducting walls ip is omitted, while in free-surface problems we also set 

a = Re(de i{ax+fSy)+ " ft ), (2.31) 

where a £ C is the com plex free-s urface amplitude. In (|2.30|) and (|2.31j) we have adopted 
the convention used bv lHol l|l989f) . under which Reft) =: r gives the modal growth rate 



whereas C := — Im(7)/(a 2 + /3 2 ) 1 / 2 is the phase velocity. The complex phase velocity 
c := i7/(a 2 -f-/3 2 ) 1 / 2 , where Re(c) = C and Im(c)(q 2 + B 2 ) 1 / 2 — r is frequently employed 
in the literature (e.g. lYihlll963l fl969l : iTakashima 1996f l in place of 7. 



Substituting ()2.30[) into (|2.7p and (|2.5I 6) leads to the set of coupled ordinary differential 
equations 

Ku x = iap - i(aB x + BB y )b x - B z D b x - D(B x )b z + (D U)u z , (2.32 a) 

Au y = iBp - i{aB x + [3B v )b y - B z Db y - T>(B y )b z , (2.32 6) 

Au z — Dp — i(aB x + BB y )b z - B z D b z , (2.32 c) 
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= i(au x +/3u y ) + Du z , (2.32 d) 

A m b x = -i(aB x + (3B y )u x - B z D u x + (D B x )u z - (D 17)6,, (2.32 e) 

A m b y = -i(aB x + /3B v )u v -B z Bu y + {B B y )u Zl (2.32/) 

A m b z = -i(aB x + (3B y )u z - B z Du z , (2.32 g) 

= i(ab x +Pb y ) +D« Z , (2.32 A) 

= D 2 ^- (a 2 +/3 2 )^, (2.32 i) 

where A := (D 2 -(a 2 + /3 2 )) Re^ 1 — (j+iaU) and A m := (D 2 -(a 2 +[3 2 ))Rm~ 1 -("/+iaU). 
Here, the velocity eigenfunctions are subject to the homogeneous boundary conditions 

u x (z w ) = u y {z w ) = u z (z w ) = (2.33) 

at the no-slip boundaries, which follow from (12.171 6). Moreover, if the walls are insulating, 
(|2~TTj) leads to 

b x (z w ) + ia-ip(z w ) = 0, 63, (z™) + i0i>(z w ) = 0, 
6 z (z«,) +Dii(z w ) = 0, a6j / (^) - f3b x (z w ) = 0, 

while boundary conditions for the magnetic field eigenfunctions in conducting- wall prob- 
lems are, in accordance with (|2.13[) . 

Vb x (z w ) =T>b v (z w ) = b z (z w ) = 0. (2.35) 

At the free surface, the kinematic and stress conditions, respectively (|2.18p and ()2.24|) . 
yield 

u 2 (0)-(7 + ia[/(0))a = 0, (2.36 a) 

Dii x (0) + iau z (Q) +D 2 [/(0)a = 0, Bu v (0) + i(3u z (0) = 0, (2.36 b, c) 

= p(0) - 2Re- 1 D u,(0) - (fl B (0)6 a (0) + B y (Q)&„(0) + B, (0)6,(0)) 

cos9 a 2 +f3 2 „,^ nn ,rt 2ia V ( 2 ' 36 rf ) 



TV 2 " + + B x (0)DB x (0)+B y (0)DB y (0)-—DU(0))a, 

while the insulating boundary conditions (|2.16[) become 

b x (0)+ia^(0)+DB x (0)d = 0, b y (0) + i/3$(0) + DB y (0)d = 0, (2.37 a, 6) 

6 2 (0) + D^(0) = 0, ab y (0) - (3b x (0) + {aDB y (0) - (3DB x (0))a = 0. (2.37 c, d) 

Equations (|2.32[) . in conjunction with the prescribed boundary conditions (in chan- 
nel problems these are (|2.33|) . and either (|2.34[) or (|2.35p . while in free-surface problems 
the boundary conditions are (|2~33"|1 . l[2T55]) . (|237|) . and either (|2~33|) or (|2~35|l ). consti- 
tute a differential eigenvalue problem, which must be solved for the eigenvalue 7, the 
eigenfunctions Ui, p, bi and, where appropriate, ip and/or d. As with several ot her hy- 
drody namic stability problems, it is possible to derive a Squire transformation ( Squire! 
Il933f) . mapping each three-dimensional normal mode to a two-dimensional one, with 
B y = j3 = u y = b y — 0, and smaller or equal growth rate Re(7). In the free-surface MHD 
flows studied here the Squire-transformed variables are 

a := (a 2 + 2 ) 1 ' 2 , 7 := aj/a, (2.38 a) 

Re;=aRe/a, Bm = Pm, Ga ;= Ga, Ca := Ca, (2.38 6) 

U := U, B x := (aB x + 0B v )/a, B z := aB z /a, (2.38 c) 



Large-wavelength instabilities in free-surface Hartmann flow 



11 



u x := (au x +/3u v )/d, u z := u z , b x ■= {ab x + j3b y )/a, b z = b z , (2.38a) 

p:=dp/a, a:=ad/a, ip := ip, (2.386) 

where, for reasons that will become clear below, we have opted to express the transfor- 
mations for Rm, Fr, and We implicitly through the corresponding ones for f m, Ga, and 



mations tor Km, rr, ana We implicitly tnrougn tne corresponding ones lor rm, Ga, an 
Ga. It is well known (e.g. Stuart fl954 ; Lock 1955 : Betchov fc Criminalelll967 ) that U{, b t , 



p, and tp satisfy (|2.32l a. c-e, g-i) with Hi i— > Hi, bi i— > bi, p i— » p, ip i— > V> ^ l— * ^ > ^» l— * 
a i— > a, 7 t— > 7, i?e i— ► i?e, i2?7i i— > i?m := PmRe and, importantly, /3 i— > 0. Here we verify 
that the transformation is also compatible with the non-trivial boundary conditions in 
free-surface MHD, but only if U meets the shear- free boundary condition Q2.23I 6). Specif- 
ically, making suitable linear combinations of ()2.36j) . and using (|2.37l ri). it is possible to 
derive the relations 

u z (0)-(7 + idJ7(Q))a = 0, Du x (0) + iau z (Q) + D 2 U(0)a = 0, (2.39 a, b) 

= p(Q) - 27?;- 1 D« 2 (0) - (5,(0)6,(0) + (0)6,(0)) 



(cos{0)Fr 



-2 i ~2 



<5 2 We" 1 + B,(0) D B x (0) - 2i(i?ea)^ 1 D [7(0) a, 



(2.39 c) 



where Fr := ReGa 1 I 2 and VFe := Ca/Re, while linear combinations of (|2.37|) lead to 

6,(0) + ici^(0) +DB x (0)& = 0, 6 Z (0) + D -0(0) = 0. (2.40a, 6) 

Equations (|2T551 a). (12351 6). P^DI a). and (|2~4l3l 6) are structurally similar to (|23Sl a). 
(12.361 c). (12.371 a 1 ). and (|2.37l c). respectively, and it is also straightforward to check that 
Uj and bi satisfy Squire-transformed versions of (|2.33p - (|2.35p . On the other hand, a 
comparison immediately reveals that (|2.39l c) and (|2.36l d) are compatible only if D U(0) = 
0. Thus, unlike channel problems, the validity of the transformation in free-surface flows 
depends on the functional form of the steady-state velocity. 

In plane Poiseuille flow, the correspondence between three and two-dimensional modes 
implies that for the purpose of determining the minimum (critical) Reynolds number 
for instability it suffices to restrict attention to two-djimensional normal modes. That is, 
it follows from the inequalities Re (7) ^ Re (7) and Re ^ Re, which are a consequence 
of (|2.38l a) and (|2.38l c). that to each unstable three-dimensional mode there corresponds a 
less or equally stable two-dimensional one at smaller or equal Reynolds number. However, 
in the multiple-parameter problems studied here the question of whether or not the 
critical mode is two-dimensional depends on the path followed in parameter space as 
Re is increased. In particular, if the process of increasing the Reynolds number modifies 
any of the flow parameters that are unchanged by the Squire transformation, the critical 
two-dimensional mode is not guaranteed to be of vanishi ng spanwise wavenumber. 

The latter observation has been made by iHuntJ (|l966h for channel MHD flows under 



purely streamwise applied magnetic fields (B y = B z =0), where, according to (|2.38l c). 
B x and, equivalently, the Alfven number Al = 1/B X are unchanged by the transforma- 
tion. Hence, if Al and Pm are held fixed as the Reynolds number of the three-dimensional 
problem is increased, then indeed the first mode that becomes unstable has (3 — 0. How- 
ever, if one requires the channel width, the fluid's material properties, and the external 
magnetic field strength all to remain fixed while the flow speed is increased (as was as- 
sumed by Hunt), then as Re grows Al necessarily increases as well, and three-dimensional 
modes may become unstable first. 

In the problems with flow-normal external magnetic field (B' x = B y = 0) studied here, 
(|2.38I 6) and (|2.38l c) imply that the Hartmann number Ha := Pm 1 / 2 ReB z = Ha and the 
induced magnetic field profile B := B x j ' (Pm x l 2 Ha) — B of the two-dimensional flow are 
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the same as in the three-dimensional case. Since the Hartmann number is independent of 
the characteristic velocity £/*, this in turn indicates that when channel Hartmann flow is 
driven at progressively higher speeds with all geometrical and material parameters held 
fixed (as can be accomplished by means of a pump generating the streamwise pressure 
gradient 77 (12.27[) ). the critical mode has purely streamwise wavevector. In free-surface 
problems, (|2.38I 6) necessitates that Ga and Ca are additionally constrained, but this 
cannot be accomplished simply by varying [/» with all other aspects of the problem fixed. 
This is because (i) the capillary number Ca directly depends on U* , and (ii) the boundary 
condition (|2.23l a). in conjunction with (|2.27[) . introduces the parameter interdependence 

Re = Gat&n(8)t&nh(Ha/2)/Ha or Re = Gat&n(9){scch(Ha) - 1)/Ha 2 , (2.41) 

valid respectively for problems with insulating and conducting walls. In the expressions 
above only Re depends on £7* , indicating that the steady-state velocity cannot be changed 
without modifying at least some of the remaining properties of the flow. An option 
compatible with (|2.38[) would be to fix 17* , and vary Re by changing the fluid thickness 
I, the inclination angle 6, and the external magnetic field B z in a manner that Ga and 
Ha remain constant. 

2.6. Two- dimensional normal modes 

In the two-dimensional normal-mode formulation, where j3, u y , b y , and B y are all set to 
zero, the divergence-free conditions (|2.32l ri. h) can be used to eliminate the streamwise 
velocity and magnetic field eigenfunctions, giving 

u(r, t) = Re((iD u(z)/a, 0, u{z))e iax+ ^), 6(r, t) = Re((i D b(z)/a, 0, b{z))e lax+ ^) 

(2.42) 

for the perturbed velocity and magnetic fields, where we have dropped the z subscript 
from u z and b z for notational clarity. Substituting (|2.42|) into (|2. 321 a. c, g), and eliminating 
the pressure eigenfunction then leads to the coup led OS and induction equations (e.g. 
Bctch ov fc Criminaldll967t iMuller fc Buhlerll200lh . 



i?e _1 (D 2 -a 2 ) 2 u - (7 + ia[/)(D 2 -a 2 )u + ia(D 2 U)u 

+ {iaB x + B z D)(D 2 -c?)b - ia(D 2 B x )b = 0, (2.43 a) 

and 

Rm-\D 2 -a 2 )b - (7 + iaU)b + (iaB x + B z B)u = 0, (2.43 6) 

respectively. 

Whenever applicable, we write 

b' = Re((i D b'(z) /a, 0, V ' {z))e iax+ ^) = -Re{ia^(z), 0, D iP(z)e iax+ ' 1t )). (2.44) 

for the external magnetic field perturbations. Laplace's equation (|2.32l iT) for the magnetic 
potential then becomes (D 2 —a 2 )^ = 0, which, in conjunction with the condition that rfc 
vanishes at infinity, has the closed-form solutions 

ib( z ) = {\)> ' ih{z) = {\) J ' (2.45 a. b) 

vw \^{-l)e a{z+1 \ z<-l, W \tp(-l)e^ z+1 \ z<-l, V ' 

respectively for free-surface and channel problems with insulating walls. If the walls are 
conducting, only the expression valid for z > is retained in free-surface problems, 
whereas in channel problems tf> is dropped altogether from the formulation. 

In inductionless problems, (|2.8p substituted into (|2.7l o) leads to the modified OS equa- 



Large-wavelength instabilities in free-surface Hartmann flow 13 

tion l|Stuardll954l lLocklll955h 

(D 2 -a 2 fu - (iaH x + H z D) 2 u - i?e( 7 + iaU)(B 2 -a 2 )u + iaRe(D 2 U)u = 0, (2.46) 

where H x := (ReRm) 1 / 2 B x . Equation (|2.46p , which replaces the coupled system (|2.43p , 
can also be obtained by letting Pm \ with H x and H z fixed. In that limit, the induced 
magnetic field I x vanishes (see ^2.4p , and (|2.43l M reduces to 

(iaB x + B z D)(D 2 -c?)b = -Rm(iaB x + B z D)u, (2.47) 

leading to (|2.46[) upon substitution in (|2.43l aV 

As for the boundary conditions, substituting for u x , b x , and p in the no-slip, kinematic, 
and stress conditions (|2.33[) and (|2.36[) by means of (|2.32l c. d, h), leads to 

ii(z m )=D!i(z,„) = 0, (2.48 a) 

u(0)-(7 + iaJ7(0))a = 0, D 2 u(0) + a 2 u(0) - ia D 2 U(0)d = 0, (2.48 b, c) 

and 

(((D 2 -3a 2 )D-i?e(7 + iaC/)D+iai?e(DC/))i2)| 2=0 + J Re( J B 2 (D 2 -a 2 )-ia(D J B 2; ))6| z= o 
- a 2 (GaRe^ 1 +a 2 Ca~ 1 + ReB x (0)D B x (0) - 2iaD[/(0)) a = 0. (2.48 d) 

Moreover, using (|2.45|) to eliminate -0, the magnetic field boundary conditions at insu- 
lating boundaries, ([2~34]) and ([2~37| . yield 

D 6(±1) ± ab(±l) = (2.49) 

and 

D6(-l)-a&(-l) = J D&(0)+a&(0)-iQ!DS x (0)a = J (2.50 a, 6) 

respectively for channel and free-surface problems. If, on the other hand, the walls are 
conducting (|2.35[) leads to 

b(z w ) = 0. (2.51) 

In inductionless problems, the boundary conditions for b are not required, and (I2.48l d) 
becomes 

((D 2 -3a 2 ) D -i?e( 7 + iaU) D +iaRe(D U) - H z {\aH x + H z D))i2| z=0 

- a 2 ( GaRe^ 1 + a 2 Ca~ x - 2iaDf7(0)) a = 0. (2.52) 

To summarise, in both of the free-surface and channel MHD stability problems studied 
here, the steady-state configuration and the governing differential equations are respec- 
tively (|2.26p (with z restricted to the appropriate interval, and the integration constants 
set according to the wall conductivity) and (|2.43|) . The boundary conditions for free- 
surface problems with insulating wall are (|2.48|) and (|2.50p , while if the wall is conduct- 
ing (|2.5ip is enforced in place of (|2.50l a) . In channel problems the boundary conditions 
are (|2.48l a). and either (|2.49p or (|2.5ip . In inductionless problems, the governing equa- 
tions are replaced by (|2.46p . the magnetic field boundary conditions (|2.49p - (|2.5ip are 
omitted, and (|2.48l rf) is replaced by (|2.52[) . 



3. Energy balance 

3.1. Formulation for two-dimensional perturbations 

Following the analysis bv lStuartl (|l954ft and lSmith fc Davisl (|l982l ). respectively for chan- 
nel flows with homogeneous boundary conditions and non-MHD free-surface flows, we 
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now derive energy-balance relations for normal modes in free-surface MHD with insulat- 
ing boundary conditions. The resulting formulation will contribute towards a physical 
interpretation of th e results presented in $H and can also provid e consistency checks for 



numerical schemes (S mith fc Davis! 1982 ; Giannakis et al. 2008 of) . In order to keep com 



plexity at a minimum, we restrict attention to two-dimensional normal modes, setting 
the spanwise wavenumber j3 equal to zero and assigning an arbitrary length L y to the 
size of the domain in the y direction. Moreover, we assume that the steady-state velocity 
and magnetic field profiles satisfy (|2.25p subject to the boundary conditions (|2.10l a). 
(|2.15l a). (|2.17l a). and (|2.23l fc). Similar energy equations apply for the other types of sta- 
bility problems studied here, but in the interests of brevity, we do not explicitly consider 
their derivation. In this section we work in Cartesian tensor notation, using m, bi, and b\ 
to respectively denote the components of u, b, and b' in the (x, y, z) coordinate system, 
and eijk to denote the Levi-Civita symbol. Summation is assumed over repeated tensorial 
indices. 

First, we define the kinetic energy density and the internal magnetic energy density of 
the perturbations as £ u := UiUi/2 and £ b := bibi/2, respectively. The integrals of these 
quantities over the unperturbed fluid domain fl := (0,L X ) x (0, Ly) x (—1,0), where 
L x := 2n/a is the modal wavelength, then yield the total kinetic and internal magnetic 
energies 



E u := / dV£ u , E b := / dV £ b , (3.1 a, b) 

J a J a 

where dV := dxdydz. Similarly, the magnetic energy density in the exterior of the fluid, 
given by £y := 6-6-/2, leads to the external magnetic energies 

E V -:= dV£ v , E v+ := dV£ v , E v :=E V -+E v+i (3.2) 
J n- J n + 

where Q- := (Q,L X ) x (0,L y ) x (-oo,-l) and i?+ := (Q,L X ) x [0,L y ) x (0,oo). The 
two forms of energy associated with the free surface are the potential energy E p and 
the surface-tension energy E a , defined in terms of the corresponding densities (per unit 
surface) as 



E p := / dS£ p , E a := / dS£ a , (3.3) 

JdO s JdO s 

where dft s :— (0,L X ) x (0, L y ), dS dxdy, and 

£ p := (cos(6)Fr- 2 + B x (0)D B x (0)) a 2 /2, £ a := {d x af/{2 We). (3.4) 

We remark that E p is equal to the work needed to displace the free surface from z = to 
z = a{x,y,t) in the presence of the gravitational field and the flow- normal gradient of the 
steady-state magnetic pressure. Also, noting that (1 + (c^a) 2 ) 1 / 2 dx dy = (1 + (d x a) 2 + 
0(a 4 )) dx dy is the area element on the free surface, £„ is equal to the work done against 
surface tension in order to increase the free-surface area from its unperturbed value to 
that corresponding to the amplitude a(x, y, t). The potential and surface-tension energies 
make up the total free-surface energy 

E a :=E p + E a . (3.5) 

The time-evolution of the energy in the fluid domain follows from the linearised equa- 
tions (|2.7[) . which, upon elimination of the Lorentz force / in (|2.7l a) by means of (|2.7l b). 
read 

dtUi = —Uj djUi — uj djUi + bj djBi + Re~ 1 Aiii + Bj djbi — dip, (3.6 a) 
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ji = Rm~~ eijkdjbk = e, + e ljk Ujb k + ei jk UjB k , (3.6 b) 

d t bi = —Uj djbi + bj djUi — Uj djBi + Rm^ 1 Abi + Bj djUi, (3.6 c) 

diUi = 0, dibi = 0. (3.6 d, e) 

In particular, forming the contraction of (|3.6l a) with u iy the contraction of (|3.6l fo) with 
bi, and adding the results together, leads to the energy equation 

d t (£ u + £ b ) = g R + g M + gj + g» + g n + d^-qf + 9l (em) + <?. t (mec,l) ) , (3.7) 

where 

g R := -UtUjdjUi, g M ■= hbjd-jUi, g.j := UiU^d-jBi - diBj), (3.8 a-c) 

9u ■= -SijSij/(2Re), g v := -Rmjij u (3.8 d, e) 

with Sij := d{Uj + djUi, and 

qf ~ Ui(£ u + £ b ), qf m) := e ijk bj (e k + e klm Ui b m ), (3.9 a, b) 

q^ mech ' :=uj(—(p-pb)6ij + Re~ l Sij), p b := Bibi. (3.9 c, d) 

We remark that in deriving (|3.7p we have used the divergence-free conditions (|2.7l e. f) 
(as well as the corresponding ones for Ui and Bi) and the relations 

UiAui = di(u r Sij) - SijSij/2, Rm~ 1 b i Ab i = di(€i jk bjj k ) - Rmjiji. (3.10) 

In a similar manner, Faraday's law dtb[ — ~^ijkdje' kl governing b[ and the external 
electric field e' iy in conjunction with the curl-free property £ij k djb' k — (this holds due to 
the insulating nature of the medium surrounding the fluid), yields the energy equation 

d t £y = eijkdiib'je'k) (3.11) 

for the external magnetic energy density. As for the surface energies (|3.4I) . multiplying 
the kinematic boundary condition (|2.18| by the free-surface amplitude a and suitable 
constants leads to the rate equations 

d t £ g = -U(0) d x £ g + (cos(9)Fr- 2 + B x {0) D B x (0)) au z \ z=0 , (3.12 a) 

d t £ a = -U(0)d x £ a + (d x (u z \ z=Q d x a) - u z \ z=0 d 2 x a)/We. (3.12 6) 



Each of the terms on the right-hand side of (13. 7j) has a physical interpretation. First, 
the source terms g R and gM in are respectively the Reynolds and Maxwell stresses, 
i.e. the energy transfer rate between the basic flow and the velocity and magnetic field 
perturbations. Physically, the Reynolds stress is an outcome of the mechanical exchange 
of energy occurring as the velocity perturbations transport mass within the non-uniform 
steady-state velocity field. In contrast, the energy transfer mechanism corresponding to 
gM is electromagnetic in nature. Its origin lies in the stretching/shrinking of the perturbed 
magnetic field b by the basic flow. Also, noting that djBi — diBj = Rmeji k J k , where J k := 
Rm~ x e k i rn diB m is the steady-state current, gj is interpreted as the energy transferred 
from the basic current to the perturbations, the so-called current interaction, while g v and 
g n are the viscous and resistive dissipation terms. Among the fluxes {qf\q { r\qi neCh) }i 
qf is the perturbation energy transported by the basic flow, qf m ^ is the electromagnetic 
energy flux, evaluated in the rest frame of the unperturbed fluid (recall that in the non- 
relativistic limit the electric-field perturbation in the rest frame of the fluid is e + U X b), 
and q( mechS > [ s momentum flux due to mechanical stresses acting on the fluid. 

We derive a global version of the above energy equations by integrating (|3.7p over (1, 
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and making use of the divergence theorem to reduce the volume integrals of qf\ qf , 
and q( mech *> to surface integrals. To begin, a consequence of the assumed periodicity in 
x is that j n &V diqf^ vanishes, and that the surface energies (j3. 3(1 obey 

d t E p = [ dS (cos{8)Fr- 2 + B x (0)DB x (0)) au z \ z=0 , (3.13a) 

JdOs 

d t E a = - I dS We~ l dlau z \ z= Q. (3.13 b) 

Also, integrating (|3. 1 1 [) over J?_ UJ?+, and using the magnetic field boundary conditions 
([2TT|l and (|2~T6)) . leads to 

d t E v = - [ dVd iq <f m) -D5 X (0) / dSae y \ z=0 - (7(0) / dS (b x b z )\ z=0 , (3.14) 
Jn Jan s Jan s 

while the relation 

d t (E p + E a ) = - [ dVd^r cK) - / dSau x \ z=Q (3.15) 

Jn He Jan B 

follows from the stress boundary conditions (|2.24[) . Then, integrating (|3.7p . and elimi- 
nating q\ em ^ and j mec ' 1 ) by means of (|3.14[) and (|3 . 1 5|) . we arrive at the conservation 
equation 

dtE = Gr + Gm + G.j + G v + G n + G av + G a j, (3.16) 
for the total energy E :— E u + E b + Ey + E p + E a . Here the volume terms 

G R := / dVg R , G M := / dVg M , G., := / dV g.j, 

Jn Jn Jn (31?) 

G v := / dVg u , G n := / dV g v 
Jn Jn 

are respectively the total Reynolds stress, Maxwell stress, current interaction, viscous 
dissipation, and resistive dissipation. Moreover, the surface terms 

/ dSau x \ z=0 , (3.18 a) 

Jan s 

G aJ := -RmJ y (fi) [ dS a(j v - (B x u z - B z u x ))\ z=0 (3.18 6) 

Jan B 

represent the energy transferred to the free surface by viscous and electromagnetic forces, 
respectively. In particular, noting that B x u z ~B z u x is the current induced by the velocity- 
field perturbations within the steady-state magnetic field, the quantity j y — (B x u z — B z u x ) 
in (|3.18I 6) can be interpreted as the current in duced by the steady-state fluid motion 
within the perturbed mag netic field (|Hundll966h . 

3.2. Energy balance for two-dimensional normal modes 

The results of the preceding section can be applied to the special case of the two- 
dimensional normal mode solutions. First, the expressions 

/o /-O 
dzE u (z), E b = Le 2rt J dzE b (z), (3.19) 



D 17(0) 
GaU ■ Re— 
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follow by substituting for m and 6, in (|3.1j) using (|2.42|) . where L := L x L y /Aa 2 is a 
constant, and 

E u := |Du| 2 + a 2 |u| 2 , £ b := | D 6| 2 + a 2 |6| 2 (3.20) 

are the modal kinetic and magnetic energy densities, averaged over the streamwise and 
spanwise directions. The external magnetic energy E' b (|3.2[) can be computed in terms 
of the internal magnetic field using the solution (|2.45l a) for ip to evaluate the integral 
over z, and the insulating boundary conditions (|2.34[) and (I2.37|) to substitute for the 
magnetic potential at the fluid domain boundaries. Specifically, we have 

E' b = Le 2rt Q 1 dz + J™ dz^j (| Do'(z)| 2 + a 2 \b'(z)\ 2 ) = Le 2rt a(\b(~l)\ 2 + \b(Q)\ 2 ). 

(3.21) 

In addition, the normal-mode solution (|2.3ip (with ,9 = 0) for the free-surface displace- 
ment in conjunction with (|3.3[) yields 

E p = Le 2rt a 2 (cos(9)Fr~ 2 + B x (0)D B x (0)) |d| 2 , E a = Le 2Ft : a 4 We _1 |d| 2 ! . (3.22) 

We treat the source terms on the right-hand side of (|3.16|) in a similar manner. The 
volume terms (|3.17[) become 

/o M 
dzg R (z), G M =Le 2rt J dzg M (z), 

/0 ,-0 ,-0 

dzgj(z), G„=Le 2rt J dzg u (z), G v = Le 2Ft J dzg n (z), 

(3.23) 

where 

g R :=2a(BU)Im(u*Du), g M := -2a(D U)Im(b* D 6), (3.24 a, 6) 

gj := 2aDBJm(ii , D6-b*Du), (3.24c) 

g v := -2iie _1 (| D 2 ?2| 2 - 2a 2 Re(u* D 2 u) + a 4 \u\ 2 ), (3.24 d) 

g n := -2Rm~ l {\ D 2 b\ 2 - 2a 2 Re(6* D 2 b) + a 4 \b\ 2 ). (3.24 e) 

Moreover, the surface terms (|3 . 1 8[) evaluate to 

G aU = -Le 2rt 2aRe- 1 B 2 U(0)lia(dI)u*(0)), (3.25 a) 

G aJ = Le 2rt 2aJ y (0){lm((D 2 6(0) - a 2 6(0))d*) (3.25 6) 

- 2aB z D B x (0)lm(d D u*(0)) + 2a 2 B x {Q) D £ x (0)Re(du*(0))), (3.25a) 

Finally, noting that the energy growth rate of the normal-mode solutions (I2.42j) and (|2.31[) 
satisfies d t E = 2FE, and inserting (l3~24l and ([3~25]) into ([3~T6]) . we obtain 

r = r R + r M + r, + r„ + r v + r aU + r aJ , (3.26) 

where Z].j := G[.]/2E. 

Equation (|3.26p expresses the modal growth rate as a sum of contributions from the 
various MHD energy generation and dissipation mechanisms. In the ensuing discussion, 
we oftentimes aggregate the terms on its right-hand side, writing r = E mec f l + F em , where 



Emech '■— Er + E a u + Ej + E^, E em := Em + E a j + E n (3.27 a, b) 

represent the net mechanical and electromagnetic contributions to F, respectively. More- 
over, we introduce normalised versions of the energy transfer densities (|3.24[) . writing 
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r^{z) :— L§[.j(z)/ (2E)\t—o for each of the source terms in (|3.24|) . so that j x dz /*[.] (z) = 
/].] . We also note that he corresponding energy transfer decomposition for inductionlcss 
problems can be obtained from (|3.26p by formally setting the magnetic field energies, Eb 
and E' b , the Maxwell stress as well as all terms dependent on the induced magnetic 
field to zero and, as follows from (|2.47[) . replacing (|3.24l e) with 

g v = -2i?e" 1 (i/f|Du| 2 + 2aH x H z lm(u* D u) + a 2 H 2 \u\ 2 ). (3.28) 



4. Results and discussion 



We now investigate the stability properties of the models established in restricting 
attention to problems with purely flow-normal external magnetic field (i.e. H x = and 
H z = Ha) and magnetic Prandtl number no greater than 10 -4 (including the induction- 
less limit Pm \ 0). Albeit small, the chosen upper boundary for Pm encompasses all 
known laboratory and industrial fluids. Our focus will be on the stability of travelling 
gravity and Alfven waves, neither of which are present in channel Hartmann flow. In 
addition, the behaviour of the hard instability mode, which is th e free-surface analogue 
of the even unstable mode in channel problems (|Takashimalll996l ). will be examined. All 
numerical work was carried out using a spe ctral Galerkin method fo r the coupled OS and 
induction equations for free-surface MHD (jGiannakis et all 12008 ah . 

Following a review of free surface-Poiseuille flow in §4.1( we consider in §4.21 induction- 
less problems, and then, in £ 14.31 flows at nonzero Pm. In the interest of commonality with 
the literature for channel Hartmann flow, we frequently use the complex phase velocity 
c := vy/a = C + ir/a, where C and r are respectively the modal phase velocity and 
growth rate, in place of the complex growth rate 7. The complex phase velocity will also 
be employed whenever reference is made to neutral-stability curves in the (Re, a) plane. 
In particular, we consider these curves to be the loci lra(c(Re, a)) = 0; a definition that 
does not necessarily agree on the a = axis with the equivalent one, Re(7(i2e, a)) = 0, 
in terms of 7 (see N4.3. 1|) . We denote throughout the critical Reynolds number for the 
onset of instability and the wavenumber and phase velocity of the critical mode by Re c , 
a c , and C c , respectively. 

Motivated by the discussion of the Squire transformation in fc|2.5l we have opted to 
perform our analysis using the parameter set {a, Re, Pm, Ha, Ga, Ca}, rather than, say, 
parameterising gravity by means of the Froude number Fr and inclination angle 6, and 
surface tension by means of the Weber number We. In particular, all ensuing calcula- 
tions where Re is varied (namely eigenvalue contours in the (Re, a) plane, and critical- 
Reynolds-number calculations) are performed at constant {Pm,Ha, Ga, Ca}, as under 
that condition the onset of instability is governed by two-dimensional ((3 = u y = b y = 0) 
modes. Computing eigenvalue contours at constant Ha and, where applicable, Pm is 



also the conventional choice in the literature for channel Hartmann flow (e.g. lLockl 
19551 Potter fc Kutchev 1973 ; Takashima 19961 ). However, in free-surface Poiseuille flow 
these calculations are typ i cally performed at fixed inclination angle lYihlll963l Eml 
19671 lYihlll969l : be Bruinl [1974 : lLam fc Bavazitoglu1ll986l : iFlorvan et al.lll987t ). rather 
than fixed Ga. Yet, as follows from ()2.41|) . as long as ReHaf tanh(i?a/2) <^ Ga or 
ReHa 2 /(sech(Ha) — 1) <C Ga, respectively for insulating and conducting lower wall, 8 
remains small, which is the case for several of the calculations presented here. Among 
the various dimensionless groups associated with surface tensi on, our ca pillary number 
Ga is eq uivalent to the parame t er S' = 3/(2Ca) emp l oyed by I Yihy 19631). wher e as, for 
instance, ISmith fc Davis] (|l982f ). lLam fc Bavazitoglul (|l986h . and IFlorvan et all (|l987l ) 
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Table 1 . Complex phase velocity c and free-surface energy E a (normalised by the total modal 
energy E = E u + E a ) of the 10 least stable modes of the spectrum in figure [3] The modes are 
tabulated in order of decreasing Im(c), and labelled Ai, Pi, Si, or F according to their family, 
where i denotes the rank, again in order of decreasing Im(c), within a given family. 



use S = Re/Ca, S t = 1/We = l/(ReCa), and C = 3 1/3 (Ga/ cos(6>)) 2/3 sin(0)/(2 Ca), 
respectively. 

Unless otherwise stated, we set Ga — 8.3 x 10 7 , which, for a typical liquid metal with 
kinematic viscosity v = 3 x 10~ 7 m 2 s _1 in a g = 9.81 m s~ 2 gravitational field, cor- 
responds to I 3 cos(9) = (0.00913 m) 3 . Capillary effects are of minor importance for the 
instabilities we wish to explore, but we nominally work at Ca = 0.07, which is the capil- 
lary number computed for dynamic viscosity \i — 1.5 x 10™ 3 N s m~ 2 and surface-tension 
coefficient a = 0.1 N m (both of which are typical liquid- metal values), assuming a 
velocity scale [/» = 4.7 m s -1 . 

4.1. Free-surface Poiseuille flow 

Our baseline scenario is non-MHD free-surface flow with the Poiseuille velocity profile 
U{z) = 1 — z 2 . As shown in the spectrum in figure[3]and table[TJ evaluated at Re = 7x 10 5 
and a — 2 x 10~ 3 , the eigenvalues form the characte ristic three-branch structure in the 
complex plane encountered in pla ne Poiseuille flow (jMackl [l976l ; Dongarra et al. 1996 



lKirchnerll2000t nVtelenk et al\ \2000). with the difference that, because they lack reflection 
symmetry with respect to z, the modes do no t arise as symmetric-antisymmetric pairs. 
Following standard nomenclature (jMackl ~976l ). we label the branches A, P, and S, where 



< Re(c) < (U) = 2/3 and (U) < Re(c) < 1, respectively for mode s in the A and P 



famil ies, while S modes have C = (U), asymptotically as Re(c) — > — oo (jGrosch fc Salwen 



1964h . Among the A, P, and S modes, only the ones at the top end of the spectrum carry 



appreciable surface energy. For instance, in table [U E a /E drops from 0.00652 for mode 
Pi to O(10- 7 ) for mode S 2 . 

In addition to the above shear modes, the depicted free-surface spectrum contains an 
unstable surface mode, denoted by F, which propagates downstream with phase velocity 
greater than the steady-state velocity at the free surface (i.e. Re(c) > 1). This s o-called 
soft instability is driven by viscous stresses acting on the free surface (|Yihlll963h . In the 



notation of [J3j the corresponding energy transfer rate r a u > \r„+rn\ to mode F exceeds 
the n et rate of energy dissipated by viscous and Reynolds stresses (see also Smith fc Davisl 



19821) . resulting in a positive growth rate r = r a u + T u + Fr. 

Besides mode F, whenever the speed of propagation of surface waves in the absence 
of a basic flow is large compared to the steady-state velocity, the spectr um also contains 



an upstream-propagating (Re(c) < 0) surface mode (e.g. figure 4 in iGiannakis et 
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Figure 3. Eigenvalues for free-surface Poiseuille flow at Re = 7 x 10 5 , a = 2 x 10~ 3 , 
Ga = 8.3 x 10 7 , and Ca = 0.07, showing the A, P, and S branches on the complex-c 
plane. Mode F, represented by a boldface marker, is unstable, and has growth-rate 
r — alm(c) = 3.7447 x 10~ 5 and energy transfer rates (the terms on the right-hand side 
of (pT26|) with r M , r,j, r v , and r a j set to zero) r R = -7.7875 x 10" 6 , F v = -1.1048 x 10 -4 , 
and r aU = 1.5572 x 10~ 4 . 



2008a). As Re grows, that mode joins the A branch, and eventuall y becomes unstable. 



The l atter instability, ofte ntimes referred to as the hard instability (|Linlll967t iDe Bruinl 
1974 ; Florvan et al. 1987|). is the free-surface analogue of the Tollmien-Schlichting wave 



in plane Poiseuille flow (|Linlll944j ). i.e. it is caused by positive Reynolds stress associated 
with a critical layer that develops for suitable values of the mode's phase velocity. 

As shown in figure[4](a), the growth-rate contours of the hard mode in the {Re, a) plane 
are qualitatively similar to those of the unstable mode in channel flow (IShenlll954h . In fact, 
if gravitational and surface-tension forces are decreased to zero the critical parameters 
of the hard mode appro ach the (Re r ., a r ., C c ) = (5772.2,1.021,0.264) values computed 
for plane Poiseuille flow (|Orszaglll97lh . revealing their common nature. Increasing Ga 
results in the instability region extending to progressively larger wavenumbers (see e.g. 
the Ha — results in table [5]), where the upper and lower branches of the neutral- 
stability curve Im(c) = intersect in a cusp-like manner. However, the hard mode's 
critical Reynolds number d oes not vary monotonically with th e strength of the flow- 
normal gravitational force ( De Brum] 19741 iFlorvan et aZ.I [l987 ) . In particular, for the 
Ca = 0.07 capillary number used here, Re c decreases with Ga < 10 7 (e.g. in table [2 Re c 
drops to 3711.3 for Ga = 8.3 x 10 6 ), but becomes an increasing function of the Galilei 
number for sufficiently strong gravitational fields, eventually exceeding the corresponding 
critical Reynolds number for channel flow. For instance, in figure [4](a) the hard mode's 
critical parameters are (Re c ,a c ,C c ) — (7361.0,2.815,0.184). 

As for the soft mode, it is evident from the structure of the Im(c) = contour in 
figure[2](a), which runs parallel to the log(o») axis for that its region of instability 

in the (Re, a) plane extends to arbitrarily small wavenumbers. This makes the soft mode 
amenable to study using regular perturbation theory for large wavelengths (a \ 0), 
when, in contrast, an analytic treatment of th e hard mode would r equire the full ma- 
chinery of singular asymptotic expansions (e.g. Drazin fc Reidl 2004 ) . In particular. lYihl 
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(1963) ias established that the lower branch of the soft mode's neutral-stability curve 



is the a = axis and further, that its upper branch, shown in figure [4](a), emanates 
from a bifurcation point located at (-Reb, 0) = ((bGa/8) 1 ^ 2 , 0), with corresponding phase 
velocity Cb = 2 (see also appendix lA.ip . 

As a check on the proximity of the bifurcation point (Reb-, 0) to the critical point 
(Re c , a c ) of the soft mode for the problem in figureH](a), which is strongly suggested by 
the direction of the Im(c) = contour for a <C 1, we have numerically computed the 
minimum Reynolds number Re m for instability at fixed a = 10~ 5 . The Re m — 7202.4298 
numerical result is very close to the analytically determined value Reb — 7202.4301 for 
the Reynolds number at the bifurcation point, as was also observed in a number of 
calculations with Ga £ [10 3 , 10 9 ]. Still, for our reference problem with (Ga, Ca) = (8.3 x 
10 7 ,0.07), Re m is smaller than the corresponding Reb by an amount of order O(10 -8 ), 
which in all likelihood is not due to numerics (e.g. the discrepancy did not disappear 
by increasing the polynomial degree of the discretisation scheme). This observation is 
consistent with a fourth-order asymptotic result that for any capillary number Ca there 
exists a lower bound in Ga above which d Re /da is negative on the Im(c) = contour, in 
the neighbourhood of the bifurcation point ijGiannakis et aLll2008M ). For Ca = 0.07 that 
lower bound amounts to Ga w 3.13 x 10 5 , indicating that for the problem in figure HI (a), 
a c = is not an exact statement. However, we expect the smallness of a c to render any 
unstable modes with Re < Reb irrelevant in a laboratory context, even if the true critical 
Reynolds number were to deviate significantly from Reb- 

4.2. Inductionless free-surface Hartmann flow 

In inductionless Hartmann flow, the external magnetic field influences normal-mode sta- 
bility on one hand by modifying the steady-state velocity profile (|2.26l o), therefore al- 
tering the energy transfer to the perturbations mediated by viscosity (represented by 
the terms and r a u in (|3.26|1 ) and, on the other hand, by means of the Lorentz force 
(|2.7I 6). whose only nonzero component f x :— —Ha 2 Re~ 1 u x is streamwise. The latter has 
a direct dissipative effect associated with resistivity (the energy transfer term 7^), but 
may also indirectly affect Ir, r a u, and the viscous-dissipation rate F„ by modifying the 
perturbed vel ocity field. In free-surface problems, as is the case with their fixed-boundary 
counterparts (|Locklll955t iPotter fc Kutchevlll973l : lTakashimalll996T ). the combined out- 
come of the external magnetic field is to suppress instabilities. In fact, even moderate 
Hartmann numbers (Ha ~ 3) are sufficient to shift the onset of the soft and hard insta- 
bilities to Reynolds numbers significantly higher than in non-MHD flows, in the manner 
illustrated by the eigenvalue contour plots in figure H](&). 

The critical Reynolds number, wavenumber, and phase velocity of the hard mode, 
computed numerically in figure [5] as a function of the Hartmann number Ha £ [0.1,200] 
for representative values of the Galilei number Ga/ (8.3 x 10 7 ) £ {0.1, 1, 10}, exhibit two 
distinct types of behaviour, depending on the relative strength of the gravitational and 
Lorentz forces. The first of these occurs when the Lorentz force is weak compared to 
gravity (e.g. the Ga — 8.3 x 10 8 example in figure 3 for Ha < 2), and is characterised by 
small variation of the critical parameters with Ha. As observed in the (Re, a) plane, the 
main influence of the applied field in this regime is to reduce the wavenumber bandwidth 
of the cusp- like tip of the hard mode's instability region, with little change in the position 
(Re c , a c ) of the intersection point between the upper and lower branches of the neutral- 
stability curve. Eventually, however, the tip collapses, and a c rapidly decreases towards 
the corresponding channel-flow result. For Hartmann numbers larger than that threshold 
the behaviour of the hard mode's critical parameters changes ch aracter and, as can be 
deduced by comparing table [2] to the calculations in table 1 of Takashimal ( 19961 ). it 
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Figure 4. Imaginary (left-hand panels) and real (right-hand panels) parts of the complex phase 
velocity c of the least stable mode in the (Re, a) plane, computed at constant Galilei and capillary 
numbers Ga = 8.3 x 10 7 and Ca = 0.07, for free-surface Poiseuille flow (a), and inductionless 
free-surface Hartmann flow with Ha — 3 (b). The regions of instability for the soft and hard 
modes, indicated in the left-hand panels, are distinguishable by the corresponding phase velocity 
Re(c), which exceeds unity in the case of the soft mode, but is less than the mean steady-state 
speed (17) for the hard mode. In panels (a) and (b), (U) is 2/3 and 0.742, respectively. 



becomes nearly identical to that of the unstable mode in inductionless channel Hartmann 
flow. In particular, the wavelength of the critical mo de becomes shorter, as expected 
from the decreasing thickness of the Hartmann layer (|Locklll9551 ). and for sufficiently 
strong fields the critical Reynolds number as a function of Ha is well describe d by the 
Re c = 48,250 Ha linear increase computed bv iLingwood k, Alboussiere ( 1999h for the 
unbounded Hartmann layer. In separate test calculations, where the Lorentz-force terms 
in (|2.46[) were set to zero but the Hartmann velocity profile was retained, we have observed 
that the critical parameters of the ha rd mode rem ain close to the results in table [H 
in agreement with the observation by lLockl (Il955l) that the principal contribution to 
the behaviour of (Re c , a c , C c ) comes from the modification of basic flow, rather than 
electromagnetic forces acting on the perturbed velocity field. 

As for the soft mode, it follows from the large-wavelength analysis in appendix IA.1I 
that when Ha is nonzero the a = axis remains part of its neutral-stability curve, and 
a bifurcation point (r?ef,,0), from which the upper part of the neutral-stability curve 
branches off, is again present in the (Re, a) plane. In particular, the position of the 
bifurcation point on the a = axis and the corresponding modal phase velocity C&, 
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Figure 5. Critical Reynolds number Re c , wavenumber a c , and phase velocity C c of the hard 
mode in inductionless free-surface Hartmann flow with Ga/ (8.3 x 10 7 ) £ {0.1, 1, 10}, and val- 
ues of the Hartmann number Ha logarithmically spaced on the interval [0.1,200]. The cap- 
illary number is Ca = 0.07 throughout. The critical parameters for channel Hartmann flow 
(Takashima 199 6f), as well as t he result Re c = 48,250 Ha for the unbounded Hartmann layer 
Lingwood & Alboussicrc 1999), are plotted in dashed lines for reference. 



respectively determined from the coefficients 72 and 71 in the perturbative expansion 
7 = 71 a + 72a 2 + O(a) 3 for the complex growth rate 7, are given by 

Re _ (8Ga) 1 / 2 smh(Ha/2)(Ha-tanh(Ha)) 1 / 2 
6b ~ (Hacoth(Ha/2) sech 3 (Ha) (2 Ha(2 + cosh(2#a)) - 3 sinh(2i7a))) 1 /2 ' ' a) 
C b = 1 + sech(iJa), (4.16) 



where Reb \ (SGa/8) 1 / 2 and Cb f 2 tend to their non-MHD values when Ha is de- 
creased to zero. Contrary to the non-MHD case examined in H4.ll however, for Hartmann 



24 



D. Giannakis, R. Rosner and P. F. Fischer 



Ha 












Ga 


= 8.3 x 10 6 















3.7113036E + 03 


1.87198E + 00 


2.49413E 


-01 


0.5 




4.5019913E + 03 


1.59323E + 00 


2.48917E 


- 01 


1 




8.5959446E + 03 


1.13443E + 00 


2.34231E 


- 01 


2 




2.8176098E + 04 


9.41259E - 01 


1.92059E 


- 01 


5 




1.6404990E + 05 


1.13450E + 00 


1.56426E 


- 01 


10 




4.3981065E + 05 


1.73916E + 00 


1.54789E 


- 01 


zu 




y.Dl/DDOO-Cj ~r UO 


O.Zo/Ol-tij ~r UU 


_L.OOU11.Cj 


m 

— Ul 


50 




2.4155501E + 06 


8.07657E + 00 


1.55029E 


- 01 


100 




4.8311017E + 06 


1.61537E + 01 


1.55030E 


- 01 


Ga 


= 8.3 x 10 7 















7.3610164E + 03 


2.81462E + 00 


1.84251E 


- 01 


0.5 




7.4343292E + 03 


2.77817E + 00 


1.85689E 


- 01 


1 




7.7154319E + 03 


2.64647E + 00 


1.89974E 


- 01 


2 




2.3929863E + 04 


1.10416E + 00 


1.91208E 


- 01 


5 




1.6378495E + 05 


1.13615E + 00 


1.56420E 


- 01 


10 




4.3979016E + 05 


1.73922E + 00 


1.54788E 


- 01 


20 




9.6176624E + 05 


3.23764E + 00 


1.55011E 


- 01 


50 




2.4155501E + 06 


8.07657E + 00 


1.55029E 


- 01 


100 




4.8311017E + 06 


1.61537E + 01 


1.55030E 


- 01 


Ga 


= 8.3 x 10 s 















1.9476764E + 04 


3.30597E + 00 


1.33088E 


- 01 


0.5 




1.9531440E + 04 


3.29373E + 00 


1.34239E 


- 01 


1 




1.9704003E + 04 


3.25532E + 00 


1.37611E 


- 01 


2 




2.0605875E + 04 


3.06079E + 00 


1.50076E 


- 01 


5 




1.6107058E + 05 


1.15339E + 00 


1.56349E 


- 01 


10 




4.3958470E + 05 


1.73985E + 00 
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- 01 
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- 01 
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- 01 
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Table 2. Critical Reynolds number Re c , wavenumber a c , and phase velocity C c of 
the hard mode in inductionless free-surface Hartmann flow, computed for Galilei number 
Ga/(8.3 x 10 7 ) £ {0.1,1,10}, capillary number Ga = 0.07, and representative values of the 
Hartmann number Ha in the interval [0, 100]. 



numbers lying in a relatively narrow band the bifurcation point becomes clearly sepa- 
rated from the critical point (Re c ,a c ). This is illustrated in figure [S] and table [31 where 
for the examined values of the Galilei number Ga/8.3 x 10 7 G {0.1,1,10} the critical 
wavenumber follows an a c oc Ha 3 / 4 increase for Ha < 2 (e.g. reaching a » 0.0044 for 
Ha w 1.9 and Ga = 8.3 x 10 7 ), before rapidly diminishing again at larger Hartmann 
numbers. As is the case with the Re^ oc Ga 1 ! 2 scaling in (14.11 a 1 ). all the Re c results 
nearly collapse to a single curve when scaled by Ga ' . The calculations also suggest 
that an a c oc GaT 1 ^ 2 scaling applies for the critical wavenumber, but this cannot be 
firmly confirmed with the presently attainable level of numerical accuracy and precision. 
As expected, whenever a c is small, the deviation of the critical Reynolds number and 
phase velocity from (|4.1[) is less significant, but even when a c is close to its maximum 
value the relative error is still acceptable. In table [3l for instance, Ret, overestimates Re c 
by approximately 20% when Ha — 2, while Cb underestimates C c by 4%. The influ- 




Figure 6. Critical Reynolds number Re c and wavenumber a c as a function of the 
Hartmann number for the soft mode in inductionless free-surface Hartmann flow with 
Ga/ (8.3 x 10 7 ) e {0.1,1,10}. The capillary number is Ca = 0.07 throughout. The Reynolds 
number at the bifurcation point Reb(Ha), given by (14.11 a). is plotted in dotted lines. The criti- 
cal phase velocity C c for Ga = 8.3 x 10 is shown in figure [TT1 The fluctuations present in the 
Ha < 0.2 and Ha > 3 results for a c are a consequence of the ill-conditioning of the calculation 
described in the caption to table [3] 



ence of a c > on the critical Reynolds number can also be observed by examining the 
Im(c) = contour in figure where Re decreases from 5.49 x 10 4 when a — 1CP 5 to 
4.74 x 10 4 w Re c when a = 0.0032 w a c . 

The agreement between the analytical results for (Reb,Cb) and the numerically com- 
puted values for (Re c ,C c ) steadily improves as a c (Ha) enters the decreasing phase 
(Ha > 2 in figure [6]). As such, (|4. ll a) can be used to deduce that at sufficiently large Ha 
the soft mode's critical Reynolds number Re c ~ (Ga/Ha) 1 ^ 2 exp(TJa) increases exponen- 
tially with the Hartmann number, and its critical phase velocity C c ~ l + 2exp(— Ha) de- 
creases exponentially towards unity. Similar small-a calculations (see (|A 28|) and (|A 29)) ) 
lead to the results that Reb also increases exponentially in (physically unrealistic) non- 
MHD problems with the Hartmann velocity profile, but only quadratically when the 
Lorentz-force terms in (|2.46[) are retained while the velocity profile keeps its non-MHD 
parabolic form. Therefore, the formation of the Hartmann velocity profile is the main 
driver of the critical-parameter behaviour of the soft instability as well, although it should 
be noted that the exponent in (|A 281 a) is smaller than the corresponding one derived 
from (|4.1l a) and, as a result, the critical Reynolds number of the full problem, including 
both Lorentz forces and the Hartmann velocity profile, outgrows that of the non-MHD 
test problem without bound. 

Turning now to the behaviour of the eigenvalues on the complex-c plane, a prominent 
feature of inductionless Hartmann flow, illustrated in figure [7] and movie 1 (see table |4] 
for corresponding numerical data) , is that as Ha increases the P branch of the spectrum 
becomes aligned with the S branch, and the eigenvalues in the A branch collapse towards 
the P-S branch intersection point. In addition, with the exception of mode F, which 
is seen to move along the Re(c)-axis towards Re(c) = 1, the eigenvalues are translated 
to smaller values of Im(c) (quadratically with Ha, as shown in figure [8]). As in non- 
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3.2644E + 01 
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Table 3. Critical Reynolds number Re c , wavenumber a c , and phase velocity C c of 
the soft mode in inductionless free-surface Hartmann flow, computed for Galilei number 
Ga/(8.3 x 10 7 ) £ {0.1,1,10}, capillary number Ga = 0.07, and representative values of the 
Hartmann number Ha in the interval [0.1, 8]. Also shown are the Reynolds number Reb at the 
bifurcation point and the corresponding phase velocity Cb, determined by (|4.ip . In order to 
illustrate the dependence of the critical parameters on Ga, the results for Re c and Reb have 
been scaled by Ga 1 / 2 , while a c has been scaled by Ga" 1 / 2 . We remark that because the Im(c) 
contours for mode F, which participates in the soft instability, are nearly parallel to the log(a) 
axis when a <C 1 (see figure \%\b), and the gradient of Re(7(7?e,a)) becomes shallow as Ha 
grows (this is a consequence of the strong- field neutrality of mode F discussed in the main text), 
critical-parameter calculations for the soft instability are significantly more poorly conditioned 
than the corresponding ones for the hard mode. As a result, the number of attained significant 
digits in the calculations is smaller than in table [2] especially so for a c . In addition, there is 
an evidence of an O(10 -4 ) systematic drift in the results for a c when the optimisation solver 
used to compute (Re c ,a c ,C c ) is restarted with initial conditions determined from the output 
of preceding iterations, indicating that, with the current computational resources, some of our 
results have not yet reached an asymptotic limit. However, we do not expect this to impart 
significant changes to the shape of the curves in figure [6] 



MHD problems, the phase velocity of the S-family modes is (asymptotically) equal to 
the average steady-state speed (12.29[) , which approaches unity as Ha grows, and Re(c) 
lies in the interval {(U), 1) for modes in the P branch. Moreover, the phase velocity of 
mode F remains greater than unity, even for strong fields (Ha ~ 10 3 ). By performing 
suitable test calculations, we have verified that the branch-alignment and Im(c)-decrease 
observed for the A, P, and S modes are independently caused by the formation of the 
Hartmann velocity profile and the Lorentz force, respectively. 

The exponential Reb(Ha) growth in (|4.1l a) is a somewhat misleading indicator for the 
magnetic field's stabilising effect on mode F, which participates in the soft instability. 
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Figure 7. Eigenvalues of inductionless free-surface Hartmann flow at Re — 7x 10 , a — 2x 10~ , 
Ga — 8.3 x 10 7 , and Ca — 7 x 10~ 2 . The Hartmann number in panels (a) and (b) is Ha — 10 
and 20, respectively. The evolution of this spectrum with Ha £ [0.1,50] is shown in movie 1, 
available with the online version of the paper. 



The reason is that, unlike the remaining modes in the spectrum (as well as all channel 
modes), whose decay rate increases quadratically with Ha as a consequence of Lorentz 
damping, mode F becomes asymptotically neutral for large magnetic field strengths. 
This behaviour is illustrated in figure [5J where the complex phase velocity, as well as 
results for the energy components and energy transfer rates, are plotted as a function 
of Ha £ [10~ 2 , 10 3 ] for the ten least stable modes of the non-MHD problem in figure [3l 
As the Hartmann number grows, the total mechanical energy transfer rate r mec h ■= 
+ r v + r a jj (the inductionless version of (|3. 271 a)) to mode F, experiences a sharp 
drop, caused by a decrease in F a u associated with the flattening of the velocity profile. 
This, in conjunction with resistive dissipation r v , which in inductionless problems is the 
only component of the electromagnetic energy transfer rate F em , suffices to stabilise the 
mode for all Ha > 6. However, instead of growing quadratically with Ha, as it does 
for the A, P, and S families of modes, the decay rate —T of mode F turns around, 
and approaches zero following a Ha~ 2 scaling. At the same time, the mode's energy 
content becomes almost entirely potential, with the kinetic energy following the power 
law E u /E oc Ha~ i . In contrast, for sufficiently large Hartmann numbers, free-surface 
oscillations become negligible for the A, P, and S modes, as manifested by their decaying 
surface energy E a . 

Figure [5] shows that as the Hartmann number is increased, the velocity eigenfunction 
u(z) corresponding to mode F evolves from a typical surface-wave-like profile at small 
Hartmann numbers to a state of nearly constant shear, characterised by uniform distri- 
bution of the kinetic energy away from the wall. Moreover, in line with the kinematic 
boundary condition (|2.48l b) with C w 1 (i.e. 7 « — ia), at large Hartmann numbers 
u(0) exhibits a 180° phase difference relative to the free-surface oscillation amplitude a, 
which, in the real representation corresponds to the streamwise velocity perturbations 
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Table 4. Complex phase velocity c and free-surface energy E a , normalised by the total energy 
E = E u + E a , of the 10 least stable modes of the inductionless problems in figure [7] The mean 
steady-state velocity (U), given by l|Z29|l . is 0.9001 {Ha = 10) and 0.9500 (Ha = 20). Due to 
the alignment of the P and S branches, there exists an ambiguity in distinguishing between the 
most stable P mode and the least stable S mode. Here we consider that the P branch comprises 
of the first six modes (in order of decreasing Im(c)) with Re(c) > (U). 



u x (0) being 90° out of phase with the free-surface oscillation amplitude a (recall that, in 
accordance with (|2.42l a). u x = Im(D u/a exp(7i + lax))). 

The observed scaling of the kinetic energy of mode F for strong magnetic fields is 
dimensionally consistent with a time-averaged equilibrium determined by the work done 
by Lorentz and gravitational stresses acting on the free surface. We approximate this bal- 
ance by setting / ~ g, where / ~ H a 2 Re^ 1 ] u x (0)\ and g ~ cos(8)Fr~ 2 \a\ = GaRe~ \a\ 
are respectively estimates of the Lorentz and gravitational forces. Because at large Ha 
the velocity-eigenfunction gradient D u is nearly constant over the inner part of the do- 
main, and since at large wavelengths the ratio (||mz||/||wz||) 2 ~ oT 2 is expected to be 
large, / - g leads to E u /E a ~ (K(0)| cos(6) / (Fr 2 \a\)) 2 ~ Ga/Ha 4 . The latter scaling 
is consistent with the E u /E results in figure [8] (note that E M /E s=s E u /E a for E u <C E a ), 
and in separate calculations we have checked that the E u /E a oc Ga scaling applies at 
fixed Ha. Aside from figure a \T\ oc Ha~ 2 strong-field behaviour for mode F was 
also recorded in trial inductionless problems with U (z) set to the Poiseuille profile, and 
is also expected on the basis of large- wavelength approximations (see (|A 291 c)). On the 
other hand, in non-MHD calculations with the Hartmann velocity profile, as well as 
in the corresponding small-a expansion (|A 281 c). the modal growth rate r tends to a 
Ha- independent negative value for Ha 3> 1. These observations, coupled with the di- 
mensional argument for E u /E a , suggest that the strong- field neutrality of mode F is the 
outcome of a balance between gravitational and Lorentz forces, and does not rely on the 
form of the steady-state velocity profile. 
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Figure 8. Complex phase velocity c, kinetic and surface energies, E u and E a (normalised by 
the total energy E = E u + E a ), and mechanical and electromagnetic energy transfer rates, 
r m ech and r em , for the ten least stable modes of inductionless free-surface Hartmann flow at 
Re = 7 x 10 5 , a = 2 x 10" 3 , Ga = 8.3 x 10 7 , and Ca = 0.07, showing the qualitatively different 
dependence of the F and A, P, S modes on the Hartmann number Ha. In the first and third-row 
panels, solid and dashed lines respectively correspond to negative and positive values. Besides 
mode F, the curves for modes Ai, Pi, P2, and S2 (the modes indexed according to their Ha = 
values; see table [TJ are indicated. 

4.3. Free-surface Hartmann flow at Pm ^ 10~ 4 

In low-Pm channel Hartmann flow with insulating boundary conditions the critical 
Reynolds number, wavenumber, and phase velocity are known to be well-approximated 
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Figure 9. Modulus and phase of the velocity eigenfunction u(z) (left; drawn with solid and 
dashed lines, respectively), and kinetic energy density E u (z) (right) of mode F for Re — 7 x 10 5 , 
a = 2 x 1CT 3 , Ga = 8.3 x 10 7 , Ca = 0.07, and Hartmann number Ha e {3,10,20,40}. The 
phase of the eigenfunction is computed relative to the free-surface oscillation amplitude a. 



by the inductionless scheme. In particular, the calculations by iTakashima ( 1996f ) have 
established that for Pm ^ 10 -4 the relative error incurred in Re c is less than 0.004, 
even at Hartmann numbers as high as 100, where Re c is of order 10 7 . In free-surface 
flow (again with an insulating wall), however, we find that while the critical parame- 
ters of the hard mode are equally insensitive to Pm <C 1 as in channel problems, the 
soft mode's behaviour differs markedly between the small-Pm and inductionless cases. 
Moreover, the boundary conditions now support a pair of travelling Alfven waves, the 
upstream-propagating of which may become unstable at sufficiently high Alfven num- 
bers. When conducting boundary conditions are enforced, the Alfven modes are removed 
from the spectrum, and the soft mode's critical Reynolds number becomes a decreasing 
function of Ha !3> 1. 

4.3.1. Properties of the least stable modes 

We illustrate the behaviour of the top end of the spectrum for problems with an 
insulating wall in figure [TTJ1 where contours of the complex phase velocity in the (Re, a) 
plane are plotted for the least stable mode at fixed Pm = 10~ 5 (a value lying in the 
upper end of the i-Vi-regime for liquid metals) and moderately small Hartmann number 
Ha £ [3,10]. It is evident from the proximity of the portions of the inductionless and 
Pm = 10 -5 neutral-stability curves corresponding to the hard mode, as well as the close 
agreement between the critical- Reynolds- number calculations in tables [2] and El that the 
influence of a small magnetic diffusivity on the hard instability is weak, with Re c being 
slightly smaller when Pm is nonzero compared to its value in the inductionless limit. In 
the case of the soft mode, however, prominent differences in the stability properties exist 
even at small Hartmann numbers. 

First, the structure of the Im(c) contours in figure [T0l suggest that the a — axis is 
no longer part of the neutral-stability curve Im(c) = 0, and this can be confirmed by 
means of large- wavelength perturbation theory. In particular, according to the discussion 
in appendix IA.21 free-surface Hartmann flow with an insulating wall supports, besides 
mode F, a second mode with vanishing complex growth rate 7 in the limit a \ 0, of 
magnetic origin. The zeroth-order degeneracy between mode F and this magnetic mode is 
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Figure 10. Contours of the complex phase velocity c in the (Re, a) plane for the least stable 
mode of free-surface Hartmann flow with an insulating wall at Pm = 1CP 5 , Ga = 8.3 x 10 7 , and 
Ca = 0.07. The Hartmann number Ha in panels (a-d) is 3, 5, 6, and 10, respectively. The curves 
labelled 9 in the panels for Im(c) are the Im(c) = contours of the corresponding inductionless 
problems. The neutral-stability curves for Pm = 10~ 5 , labelled N, are drawn in the Re(c) panels 
for reference. 
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Table 5. Critical Reynolds number Re c , wavenumber a c , and phase velocity C c of the hard mode 
in free-surface Hartmann flow with insulating boundary conditions, computed at Galilei number 
Ga = 8.3 x 10 7 , capillary number Co = 0.07, magnetic Prandtl number Pm £ {10 -5 ,10 -4 }, 
and Hartmann number Ha £ [0.5, 100]. 



broken at first order in a, where there exist two distinct solutions, respectively 7} and 
j[ M \ for the coefficient 71 in the perturbative series 7 = 70 + 71a + 72a 2 + 0(a 3 ). Both 
solutions have negative real part for all Ha > and Pm > (Re(7^ M ') is negative even 
when Ha equals zero) , which implies that for any Reynolds number there exists an upper 
bound a m in a, below which mode F is stable. That is, Re(7) is negative for < a < a m 
or, as observed in figure fTTJl Im(c) < for ^ a < a m (cf. inductionlcss flow). We remark 
that because 7 vanishes in the limit a \ 0, it is important to distinguish between the 
definitions Im(c) = and Re(7) = for the soft mode's neutral-stability curve, since 
with the latter definition the a = axis remains part of the curve even when Re(7i) ^= 0. 
In separate numerical calculations we have checked that a m becomes smaller when Pm 
is decreased at fixed Ha, which is consistent with the large- wavelength result (|A 40[) that 
in the limit Pm \ 0, 7^ reaches the inductionless value (|A 22|) . while 7^ is singular. 

Performing a similar type of analysis establishes that (i) when the induced magnetic 
field B is set to zero, Re(7i) is still negative for mode F, as well as for the magnetic 
mode, (ii) when U and B are both set to zero, 71, now given by (|A 41[) . is negative for 
the magnetic mode but vanishes for mode F, and (iii) in Hartmann flow with a conducting 
wall, the b(— 1) = constraint imposed on the magnetic field eigenfunction eliminates 
the magnetic mode, and 7} becomes equal to the corresponding expansion coefficient 
in the inductionless limit. It therefore appears that the suppression of the soft instability 
for a \ is the combined outcome of the boundary conditions, which allow for the 
presence of the magnetic mode, and the coupling between mode F and the magnetic 
mode provided by the steady-state flow. 

The soft mode's departure from inductionless behaviour is also prominent at larger val- 
ues of aRe. As can be seen in figure flOl at moderate Hartmann numbers (Ha ~ 5) regions 
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of stability emerge in the (Re, a) plane that would contain unstable modes in the induc- 
tionless limit. Moreover, as Ha grows, a wedge-like instability region forms, extending 
to Reynolds numbers significantly smaller than in the corresponding inductionless prob- 
lems. Unstable modes may now have phase velocity smaller than unity, but Re(c) > 1 
is found to apply for the modes close to the tip of the wedge-like region (including the 
critical mode), at least for the Ha ^ 10 3 interval covered in our calculations. 

The critical parameters of the soft mode as a function of Ha £ [10 _1 , 10 3 ] are plotted 
in figure [TT] for logarithmically spaced values of the magnetic Prandtl number Pm € 
{10~ 6 , 10~ 5 , 10 -4 } (see also table [6]). As expected from the contour plots in figure [TOl 
when the Hartmann number is small, Re c is close to its value in the inductionless limit. 
For instance, the relative difference between the Ha = 1 results in table [6] and the corre- 
sponding inductionless results in table [3] is less than 3.5%. However, once the magnetic 
field strength exceeds a threshold, which decreases with Pm, Re c (Ha) branches off from 
the exponential growth (|4.1| q), and follows closely the power law Re c oc Ha 2 ' 3 . During 
that transition, the decrease of the wavenumber with Ha observed in the inductionless 
limit is reversed, switching over to an a c oc Ha 5 ^ 4 scaling. Moreover, the exponential 
decrease (|4.1I 5) of the critical phase velocity relative to the steady-state flow at the free 
surface becomes a C c — 1 oc Ha~ 4 ^ 3 power law. In this intermediate Hartmann-numbcr 
regime, the results for Re c (Ha), a c (Ha), and C c (Ha) — 1 collapse to nearly single curves 
if scaled by Pm 1 / 3 , Pm" 1 ^ 3 , and Pm~ 2 / 3 , respectively (though the agreement is not as 
good for the Pm = 10~ 4 data for C c — 1). The power-law behaviour of the critical pa- 
rameters is only transient, however. Eventually, at sufficiently large Hartmann numbers, 
Re c (Ha) levels off (for Pm = 1CP 4 this occurs around Ha = 500), and a c (Ha) becomes 
a decreasing function once again. 

The deviation of the critical Reynolds number of the soft instability from (14.11 a) is 
even more pronounced in problems with a conducting wall. As outlined in appendix I A. 2 1 
in this case the first-order coefficient in the perturbation series for 7 is given by the 
same expression as in inductionless flows (i.e. (|A 2210 . which has vanishing real part, 
and as a result the a = axis remains part of the neutral-stability curve. This enables 
the derivation of the closed-form expression (|A 42|l for the Reynolds number Reb at the 
bifurcation point, which, in the manner shown in figure fTTt becomes a decreasing function 
of Ha, varying like Reb ~ (Ga/Pm^^Ha -1 at large Hartmann numbers. Even though 
we have not explicitly computed the soft mode's critical parameters for conducting-wall 
problems, we have verified in eigenvalue contour plots that, as in inductionless flows, 
Reb is close to Re c . In any case, since Reb is an upper bound for Re c , the behaviour 
of Reb(Ha) suffices to conclude that in free-surface Hartmann flow with a conducting 
wall the external magnetic field leads to a reduction of the critical Reynolds number for 
instability for all Pm > 0. 

4.3.2. The role of the steady-state induced magnetic field 

Except for the large-wavelength instability suppression observed in problems with in- 
sulating boundary conditions, where, as stated in N4.3.11 the modal decay rate is nonzero 
to linear order in a even when B vanishes, the discrepancy between inductionless and 
nonzero- behaviour for mode F is mainly caused by the steady-state induced magnetic 
field (|2.26I 6) and the associated contribution fj := RmJ x()~ Pm 1 ! 2 Ha\ \ D£?||||&|| to 
the linearised Lorentz force Q2.7I 6) from the spanwise current J := i?m _1 V X B = 
Rm~ x HaPm 1 / 2 T> By. As a demonstration, in figure [T3l and table [7] we have computed 
the complex phase velocity of mode F, as well as certain energy components and energy 
transfer rates, as a function of Ha <E [10~ 2 , 10 3 ] for (i) Pm = 10 -5 flows with insulating 
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Figure 11. Critical Reynolds number Re c , wavenumber a c , and phase velocity C c of the soft 
mode as a function of the Hartmann number Ha G [10 _1 , 10 3 ] for inductionless problems (curves 
labelled I), and insulating-wall problems with Pm £ {10 -6 , 10 -5 , 10 -4 }. Also shown, in dotted 
lines, are the Reynolds number Ret, and the corresponding phase velocity Cb at the bifurcation 
point of the neutral-stability curve, evaluated for conducting-wall problems by means of the 
large- wavelength results (IA 42[l and 1)4.11 6). The Galilei and capillary numbers are Ga = 8.3 X 10 7 
and Ca = 0.07 throughout. As in figure [6l the fluctuations in the a c (Ha) curves are due 
to the ill-conditioning of critical-parameter calculations for the soft mode. Aside from these 
fluctuations, we do not expect significant effects on the shape of the curves due to ill-conditioning. 



and conducting walls, (ii) the corresponding inductionless problems, and (iii) insulating 
and conducting-wall problems with Pm = 10~ 5 and B set to zero. 

The results for the test problems with B = agree fairly well with the corresponding 
ones in the inductionless limit. Here the main differences are that the total electromag- 
netic energy transfer rate r em := +i~M now includes a small, positive Maxwell stress, 
and for Ha > 30, that the modal phase velocity is less than the free-surface steady-state 
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Figure 12. Complex phase velocity c, kinetic and surface energies, E u and E a (normalised by 
the total energy E), and mechanical and electromagnetic energy transfer rates, r mec u and r em , 
for mode F in free-surface Hartmann flow at Re = 7 x 10 5 , a = 2 x 1CT 3 , Ga = 8.3 x 10 7 , 
Ca = 0.07, and Ha G [10 -2 , 10 3 ]. The curves labelled (a) correspond to inductionless flow, while 

(b) and (c) were evaluated at Pm — 10~ 5 , respectively with insulating and conducting boundary 
conditions. Curves (d) and (e) are for the same problems as (b) and (c), respectively, but with 
the magnetic field profile B set to zero. In the first and third-row pan els, solid and dashed line s 
respectively correspond to negative and positive values. As discussed in Giannakis et al. (2008a), 
accurate numerical evaluation of some of the energy transfer rates can be problematic, especially 
when a large polynomial degree is required for the spectral solution to (|2.43p to converge. For 
this reason, we have not been able to evaluate r mK j and r em for Hartmann numbers as large 
as 10 3 . Roundoff errors also limit the convergence of the eigenvalue for mode F to about five 
significant digits when Ha becomes large in cases (6-e). Apart from the Re(c) graph for case 

(c) , this level of numerical precision is not sufficient to continue the logarithmic plots in the first 
row to Ha > 100, where both of |Im(c)| and |Re(c) — 1| become small compared to |c ~ 1. 
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10" 



7.2341E + 03 
9.8246E + 03 
4.7810E + 05 
2.1205E + 06 
9.3869E + 06 



5.353E 
3.069E 
9.369E 
1.670E 
2.113E 



04 
03 
04 
02 
01 



9.9504E 
7.0268E 
1.2076E 
5.7707E 
1.8874E 



01 
01 
03 
05 
06 



Pm -- 
0.1 
1 

10 

100 

1000 

Pm 
0.1 
1 

10 

100 

1000 



10" 



= 10" 



7.2343E + 03 
9.8253E + 03 
2.1891E + 05 
9.7469E + 05 
3.8912E + 06 



7.2322E + 03 
9.5341E + 03 
1.0034E + 05 
4.4148E + 05 
1.3161E + 06 



5.355E 
3.072E 
1.978E 
3.552E 
2.244E 



6.005E 
2.979E 
3.571E 
7.375E 
6.202E 



04 
03 
03 
02 
01 



04 
03 
03 
02 
02 



9.9505E 
7.0314E 
4.9470E 
2.5559E 
4.1975E 



9.9563E 
7.3307E 
1.5546E 
9.8053E 
1.4324E 



01 
01 
03 
04 
06 



01 
01 
02 
04 
06 



Table 6. Critical Reynolds number Re c , wavenumber a c , and phase velocity C c of the 
soft mode in free-surface Hartmann flow with insulating boundary conditions, computed at 
Galilei number Ga = 8.3 x 10 7 , capillary number Ca = 0.07, magnetic Prandtl number 
Pm € {10 -6 , 10 -5 , 10 -4 }, and Hartmann number Ha G [10 _1 ,10 3 ]. These calculations are 
affected by a similar ill-conditioning as the corresponding ones for inductionless problems in 
table In particular, we have observed an O(10~ 4 ) systematic drift in some of the results for 
a c , arising when the optimisation solver used to compute (Re c ,a c ,C c ) is restarted using the 
output of previous calculations for initialisation. With the presently available computational 
resources we were not able to perform a sufficiently large number of iterations so as to establish 
convergence in a c . However, we expect the results for Re c and C c to be affected by this issue to 
a lesser extent. 



velocity. The error is particularly small for the conducting-wall problem, since in that 
case the energy transfer to the modes via Maxwell stress, which is not captured by the 
inductionless scheme, is suppressed due to the nature of the boundary conditions. In con- 
trast, when B is included, the current-interaction term rj causes the total mechanical 
energy transfer rate F mec h := F v + Fr + F a u + Fj to remain positive for Hartmann num- 
bers larger than the Ha w 9 value above which it becomes negative in inductionless and 



F 



aJ: 



B = problems. The total electromagnetic energy transfer rate r em := r v + 7\f 
where the surface term r a j is appreciable and positive, also deviates markedly from its 
-ffa-dependence in the inductionless limit. 

In the insulating-wall example, the combined effect produces a more than two-fold 
increase of the Hartmann number required for instability suppression relative to the in- 
ductionless flow. Moreover, the modal energy in the strong-field limit, instead of becoming 
almost exclusively gravitational, is split into a £/a-independent mix between magnetic 
and surface parts. Although we have observed a scaling of the form Eb/E a oc Pm/ Ga for 
Ha 1, we have not been able to account for it by invoking a work-balance argument 

(cf. S3. 

If now the wall is conducting, the current interaction term becomes sufficiently 
large so as to cause the growth rate Im(c)a to asymptote to a positive value, rather than 
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Table 7. Growth rate and phase velocity, velocity and magnetic field eigenfunctions at the free- 
surface, energy components, and energy transfer rates for mode F at Re — 7 x 10 5 , a — 2 x 10 -3 , 
Ga = 8.3 x 10 7 , Ca = 0.07, Ha = 20, and Pm = 10~ 5 . The data in column 1, counting from 
the left, are for the corresponding inductionless problem, and therefore all entries involving 
the magnetic field eigenfunction are omitted. The results in columns 2 and 4 are for insulating 
boundary conditions, while for those in columns 3 and 5 conducting boundary conditions have 
been imposed (i.e. the external magnetic energy component Et- is omitted). For the calculations 
in columns 4 and 5 the induced magnetic field B has been set to zero, and as a result i~j and 
r a j vanish. 



approach zero from below. In addition, the energy in the magnetic degrees of freedom 
dominates, with both surface and kinetic contributions decaying to zero. The greater 
influence of the steady-state current in problems with a conducting wall is consistent 
with the fact that |B(z)| is of order unity throughout the core of the flow domain, while 
it is of order 1/Ha when the wall is insulating (see £ )2.4p . 

The data for -Tr and T v in table [7] exhibit a particularly large discrepancy between the 
nonzero- £? problem with an insulating wall and its inductionless counterpart, signalling 
that the Lorentz force associated with J causes significant changes in structure of the 
velocity eigenfunction u{z). Indeed, as illustrated in figure Q2J the perturbed velocity 
field of the full MHD problem bears little resemblance to the inductionless examples in 
figure [5J In particular, instead of evolving towards a state of constant shear as Ha grows, 
|w(z)|, as well as the kinetic energy density E u {z), develop a maximum in the wall region, 
where \J\ is concentrated (see figure [5]). The increased eigenfunction curvature leads in 
turn to higher viscous and resistive dissipation, but these are more than counter-balanced 
by positive current interaction, and by positive Reynolds and Maxwell stresses. As for 
the magnetic field eigenfunction b(z), its modulus varies by less than 1CP 2 over the fluid 
domain, and since the modal wavenumber a — 0.002 is small, the energy of the magnetic 
field penetrating into the exterior region, Ey (|3.2ip . exceeds the internal magnetic energy 
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Eh (|3.19|) by two orders of magnitude (see table [7j . We remark that in problems with 
a conducting wall, where the current distribution over the inner part of the domain is 
nearly uniform, u(z) is qualitatively similar to the inductionless case, and b(z) varies 
nearly linearly from zero at the wall to its free-surface value. 

4.3.3. Travelling Alfven modes 

Having studied the behaviour of the least stable modes in some detail, we now examine 
the influence of a small, but nonzero, magnetic Prandtl number on the more stable 
modes. As illustrated by the spectra in figure fT4l and table computed for Ha = 0.1 
and Pm £ {10 -5 , 10~ 4 }, when the applied magnetic field is weak, free-surface flows with 
an insulating wall exhibit, in addition to the A, F, P, and S modes, a mode labelled M 
that is not present in inductionless flows. Mode M is stable, and like S-family modes, 
its phase velocity is close to the mean steady-state speed. However, it stands out from 
the remaining modes because of (i) the strong Pm-dependence of its decay rate, which 
decreases by an order of magnitude between the Pm = 10~ 5 and Pm = 10~ 4 calculations 
in table [8] (the corresponding relative variation for the A, F, P, and S modes is less than 
10 -4 ), and (ii) its mostly magnetic energy content (e.g. in table [H Eb/E w 0.95 for 
mode M, whereas Eb/E < 10 -4 for modes in the A, F, P, and S families). Although we 
have not confirmed this analytically, numerical calculations strongly suggest that mode M 
is singular in the inductionless limit Pm \ 0. In particular, as shown in movie 2, when Pm 
is increased from small values the complex phase velocity of mode M is seen to approach 
the upper part of the spectrum from arbitrarily small values of Im(c), moving along the 
S eigenvalue branch. We remark, however, that the complex phase velocity does not cross 
the Im(c) — axis. Instead, if Pm is allowed to be of order unity, mode M eventual ly 
participates in the formation of magnetic eigenvalue branches ( Giannakis et al. 12008 ah . 



The existence of modes of magnetic origin (hence the designation M) in the spectrum 
of the coupled OS and induction equations (|2.43|) is consistent with the fact that the 
limit Pm \ at fixed Ha, which, as discussed in £|2.6i leads to the approximate stability 
equation (|2.46|) . is singular (effectively, the differential order of the stability problem is 
reduced from six to four). Modes of this type are also present in channel Hartmann flows, 
as well as in test problems with B = 0. However, we have found no evidence of mode M 
in numerical calculations with conducting boundary conditions, which correlates with 
the absence of the 7} magnetic solution in the large-wavelength approximations for 
conducting- wall problems in appendix IA. 21 

As shown in the spectra in figure 1151 evaluated for insulating-wall problems with 
Pm = 10~ 5 , Ha £ {10,20,50}, and otherwise the same parameters as the inductionless 
calculations in figure UJ when the external magnetic field is of appreciable strength the 
collapse of the A eigenvalue branch and the alignment of the P and S branches observed 
in the inductionless limit (see §4.2p are also present in Pm > problems. In addition, 
according to the data in table El the magnetic energy for all but the first handful of modes 
is relatively small (Eb/E < 10 _1 ), and in those cases the accuracy of the inductionless 
approximation is very acceptable (for instance, a comparison with table [4] shows that 
the relative error for mode P4 at Ha — 10 is at the 1% level). Still, as already discussed 
in §4.3.11 and §4.3.2[ nonzero- Pm spectra deviate substantially from the corresponding 
inductionless ones in the behaviour of mode F, whose main distinguishing features in 
figure [12] and table [5] in comparison with figure UJ and table 2] are (i) sub-unity phase 
velocity C f» 0.9959 when Ha = 10, (ii) positive growth rate in both of the Ha — 10 and 
Ha = 20 calculations, and (iii) magnetic-energy predominance for strong background 
fields (e.g. E b /E rj 0.93 for Ha = 50). 

Aside from the discrepancies associated with mode F, however, the spectra in figure 1151 
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Figure 13. Velocity eigenfunction (a), magnetic field eigenfunction (b), kinetic and magnetic 
energy densities (c), current interaction (d), viscous dissipation and current-interaction densities 
(e), and resistive dissipation and Maxwell-stress densities (f) for mode F in insulating- wall 
problems with Pm — 1CF 5 , Ha G {3, 10, 20, 40}, and remaining flow-parameters equal to those 
in the inductionless problems in figure [9] In panels (a), (b), (e), and (J) solid and dotted lines 
are associated with the lower and upper horizontal axes, respectively. In order to better depict 
the peaks in E u (z), the corresponding axis has quartic scaling. Likewise, quadratic scaling has 
been used for the z coordinate in panels (c-/), as well as the axes for r v (z) and r v (z). 
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Figure 14. Eigenvalues of free-surface Hartmann flow with an insulating wall at Re — 7 x 10 5 , 
a = 2 x 10~ 3 , Ga = 8.3 x 10 7 , Ca = 0.07, and Ha = 0.1. The magnetic Prandtl number Pm in 
panels (a) and (6) is 10 -5 and 10 -4 , respectively. Mode M, plotted using a + marker, is singular 
in the inductionless limit Pm \ 0, while mode F, plotted in boldface, is unstable in both panels. 
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Table 8. Complex phase velocity c, free-surface energy E a , and magnetic energy Eb of the 10 
least stable modes of the spectra in figure IT4l The energies E a and Eb have been normalised by 
the total energy E — E u + Eb + E a . 
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Figure 15. Eigenvalues ol free-surface Hartmann flow with insulating wall at Re — 7 x 10 5 , 
a = 2 x 1CT 3 , Pm = 10~ 5 , Ga = 8.3 x 10 7 , Ca = 0.07, and Hartmann number Ha = 10 (a), 
20 (b), and 50 (c). The modes labelled L_ and L + , and plotted using * markers, are travelling 
Alfven waves, and can be continuously traced to modes A2 and Pi, respectively, as Ha \ 0. 
The mode plotted using a + marker originates from mode M in figure [T3T a). In panels (a) and 
(6), mode F, plotted in boldface, is unstable. The evolution of this spectrum with Ha £ [0.1, 50] 
is shown in movie 3, available with the online version of the paper. 



differ from those in figure [4] in that they contain two travelling modes, labelled L_ and 
L_|_, which have no counterparts in the inductionless limit. As illustrated in movie 3, these 
modes are the outcome of a mode-conversion process, whereby the weak-field modes A2 
and Pi separate from the A and P eigenvalue branches when Ha is increased from 
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Table 9. Complex phase velocity c, surface energy E a , and magnetic energy Eh of the 10 least 
stable modes of the spectra in figure 1151 The energies E a and Eb have been normalised by the 
total energy E = E u + Eb + E a . If Ha is decreased to 0.1, mode P4 (in figure [151 plotted using a 
+ marker) can be continuously traced to mode M in the Pm = 10 -5 part of table[8] The average 
steady-state velocity (U) is, in accordance with (|2.29[) . 0.9001, 0.9500, and 0.9800, respectively 
for Ha = 10, 20, and 50. Due to the alignment of the P and S branches, the classification of the 
sixth least stable mode for Ha = 10, the eighth least stable mode for Ha — 20, and the ninth 
least stable mode for Ha = 50 as members of the P family is somewhat arbitrary; these modes 
could equally be treated as members of the S branch. 



small values, and move towards nearly symmetric positions on the complex plane about 
Re(c) — 1. At the same time, mode M, joins the P eigenvalue branch, and having lost 
the majority of its magnetic energy, behaves in the same manner as the remaining P and 
S modes. A principal feature of the L modes for sufficiently strong magnetic fields, which 
can be observed in the Ha = 20 and Ha = 50 results in table [9] is a near equipartition 
between magnetic and kinetic energy. Because of this, we classify these modes as travelling 
Alfven waves. In separate numerical calculations, we have confirmed that modes of this 
type are also part of the spectrum when the background fluid is at rest. In channel 
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problems with insulating walls, the pair of travelling; w aves becomes replaced by a single 
frozen- in magnetic mode (jBetchov fc Criminald 1967L §IX), whose phase velocity and 
kinetic energy are nearly equal to unity and zero, respectively. However, as with mode M, 
both of the travelling and frozen-in Alfven modes are removed from the spectrum when 
conducting boundary conditions are enforced. 

When a steady-state flow is present, the upstream-propagating wave L_ may become 
unstable if the Alfven number Al = RePm 1 ^ 2 / Ha of the flow is sufficiently large. This 
situation is illustrated by the eigenvalue and energy calculations in figure[lj)]and tablefTOl 
evaluated at Re = 6.3 x 10 6 , a = 1(T 4 , Pm = lfr 4 , and Ha e [1CT 2 ,10 3 ], where 
the latter Hartmann-number interval corresponds to an Alfven number decrease from 
6.3 x 10 6 {Ha = 10 -2 ) to 63 (Ha = 10 3 ). The resulting evolution of the eigenvalues in 
the complex-c plane is shown in movie 4. 

When the Hartmann number is large (Ha > 200), the travelling Alfven modes are 
clearly distinguishable by the linear _ffa-dependence of their decay rate |Im(c)a| and their 
phase velocity relative to the free-surface steady-state velocity, Re(c) — 1. The nearly equal 
split between kinetic and magnetic energies remarked upon above is also evident, both in 
figure [TBI and in the Ha — 200 results in table [TJJ] The latter also indicate that, in light 
of the smallness of the wavenumber employed in the calculation, the magnetic energy 
contained in the exterior region exceeds the internal one by more than two orders of 
magnitude. Interestingly, a duality between the mechanical and electromagnetic energy 
transfer rates appears to apply between the downstream and upstream waves, in the 
sense that the mechanical energy transfer rate r mec h = —1.129 x 10~ 4 for mode L+ 
(the downstream-propagating mode) is nearly equal to the T em term for mode L_ and 
likewise, the r mec h — 9.830 x 10 -5 term for mode L_ is close to the electromagnetic 
energy transfer rate r em for mode F. In both cases, however, the energy transfer rate 
due to Maxwell stress and the current interaction, respectively Fm and Fj, as well as the 
surface term r a j, are positive, whereas the energy transfer rate due to Reynolds stress, 
Ifl, is negative. 

At smaller Hartmann numbers (equivalently, larger Alfven numbers), the upstream- 
propagating mode tends to be advected along the direction of the basic flow, and for 10 < 
Ha < 20, the positive mechanical energy transfer rate, driven in part by positive Reynolds 
stress (see the Ha = 15 calculation in table [T0|) . exceeds the rate of energy dissipated 
electromagnetically, resulting in an instability. As Ha is further decreased, mode L_ 
is stabilised and becomes converted to mode M, i.e. its energy becomes predominantly 
magnetic, and its phase velocity C approaches the mean value of the steady-state flow. 
This mode conversion is different from the one observed in figure [15] and movie 3, where 
mode L_ develops from weak- field mode A2, indicating that the precursors of the Alfven 
modes in the weak-field spectrum are not universal. 

As for the downstream-propagating wave L + , instead of joining the P eigenvalue branch 
when Ha is decreased (which would be the case for the flow parameters in figure [15]) , 
becomes converted, around Ha — 3, to mode F, i.e. it has greater than unity phase 
velocity, and for the chosen values of the Reynolds and Galilei number, is unstable for 
Ha = 0. At the same time, the strong-field (Ha > 80) mode in figurefTBl labelled F, which 
exhibits the asymptotic neutrality that we ascribed to mode F in figureE] turns into mode 
Pi. This exchange of identity between modes F and Pi takes place at Reynolds numbers 
greater than the one employed in our preceding examples at (Re, a) = (7x 10 5 , 2 x 10~ 3 ), 
and as can be seen in figure [16] and movie 4, is accompanied by a small Hartmann- 
number interval where both modes are stable. The latter gives rise to regions of stability 
in the (Re, a) plane which would contain unstable modes in the inductionless limit (see 
figure [TO]). 
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Figure 16. Complex phase velocity c, magnetic and surface energies, E\, and E a (normalised by 
the total energy E — E u + Eh + E a ), and mechanical and electromagnetic energy transfer rates, 
r m ech and -Tem, for the ten least stable modes of free-surface Hartmann flow with insulating 
boundary conditions at Re = 6.3 x 10 6 , a = 1CT 4 , Pm = 1CT 4 , Go = 8.3 x 10 7 , Ca = 0.7 and 
Ha e [1CT 2 , 10 3 ]. In the first and third -row panels, solid and dashed lines respectively correspond 
to negative and positive values. The weak-field (Ha — 0.1) modes F, M, and Pi respectively 
become converted to modes L+, L_, and F as Ha grows. The corresponding ffa-evolution of 
the eigenvalues in the complex plane is shown in movie 4, available with the online version of 
the paper. 
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Table 10. Growth rate and phase velocity, velocity and magnetic field eigenfunctions, energy 
components, and energy transfer rates for modes F, L_, L+, M, and Pi in figure [16] for repre- 
sentative values of the Hartmann number in the interval [0.02,200]. For the reasons stated in 
the caption to figure [121 the relative numerical error in the energy-transfer-rate calculations for 
mode F at Ha = 200 is of order 0.3. 



The results for mode F in figure [TBI also show that the \T\ oc Ha~ 2 strong-field scaling 
observed in inductionless problems (e.g. figured]) does not necessarily apply in nonzero- 
Pm flows. In particular, while no such evidence was observed in the calculations in 
figure [T3l for {Re,a,Pm) = (7 x 10 5 ,2 x lCr 3 ,lCT 5 ), the growth rate of mode F for 
(Re, a, Pm) = (6.3 x 10 6 , 10 -4 , 10~ 4 ) follows an inverse-cubic decrease. This is most likely 
caused by the steady-state current, since, as we have numerically confirmed, setting B 
to zero restores the \r\ oc Ha~ 2 scaling. 
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Before closing, we note that the A, P, and S modes that do not participate in the 
modal interactions described above also exhibit new aspects of behaviour compared to 
inductionless problems. In particular, as the magnetic field strength grows their phase 
velocity experiences a series of oscillations about C — 1, of diminishing amplitude (cf. 
the inductionless calculations in figure [8l where C monotonically approaches unity from 
below), and so does their magnetic energy as it settles towards a f/a-independent equi- 
librium value. An intricate pattern of oscillation is also observed for the mechanical 
energy dissipation rate -r mec ^, which now exhibits an increasing trend with Ha 3> 1 
instead of asymptoting towards i/a-independent values. However, at least for the exam- 
ined Hartmann-number interval, |-T mec h| remains small compared to the rate of electro- 
magnetic energy dissipation — -T em , which, dominated by the resistive term i"^, grows 
quadratically with Ha. As a result, the decay rate — r = —(r mec h + r em ) of the A, P, 
and S modes exhibits to a good approximation a f oc Ha 2 scaling for strong applied 
fields, as it does in the inductionless limit. 



5. Conclusions 

A numerical investigation of the stability of temporal normal modes of free-surface 
Hartmann flow at low magnetic Prandtl numbers (Pm 10~ 4 , including the induc- 
tionless limit Pm \ 0) has been presented. Our main objective has been to study the 
influence of a flow-normal magnetic field (of associated Hartmann number Ha ^ 1000) on 
the soft and hard instability modes p resent in free-surface Poiscuillc flow (|De Brumlll974l : 



Florvan et aZ.lll987l : lYihlll963l . ll969h . imposing either insulating or conducting boundary 
conditions at the wall. 

We have confirmed that the Squire transformation for MHD ( Betchov fc Criminalel 



1967T ) is compatible with the kinematic, stress, and magnetic field-continuity boundary 



conditions for free-surface problems, but found that unless the flow is driven at constant 
capillary and Galilei numbers, respectively parameterising the surface tension and the 
flow-normal gravitational force, the onset of instability as the Reynolds number Re grows 
is not necessarily governed by a two-dimensional mode. 

In inductionless flows, we have observed that the critical Reynolds number Re c of both 
of the hard and soft instability modes increases monotonically with Ha. In particular, 
except for applied fields sufficiently weak for gravity to dominate over Lorentz forces, the 
hard mode's critical Reynolds number, as well as its critical wavenumber a c and phase 
velocity C c , were found to be very close to the correspondin g parameters of the even 
unstable mode in channel Hartmann flow ( Takashima 19961 ). reflecting the common, 
critical layer, nature of these two instabilities. In fact, for sufficiently large Hartmann 
numbers, the critical Reynolds number of both is well approximated by the linear power 
law Re c (Ha) = 4 8,250 Ha computed for the crit ical Reynolds number of the unbounded 
Hartmann layer ( Lingwood fc Alboussierel "l999h . 

As for the soft mode, we derived, using large-wavelength perturbation theory, closed- 
form expressions for the Reynolds number Rei, and phase velocity Ct, at the bifurcation 
point (i?efc, 0) in the (Re, a) plane, where a is the wavenumber, from which the upper and 
lower branches of its neutral-stability curve emanate. Even though in inductionless Hart- 
mann flow the soft mode acquires a nonzero critical wavenumber, numerically computed 
values for Re c and C c were found to be in close agreement with the analytical results for 
Reb and respectively. From these it follows that Re c ~ (Ga/ Ha) 1 / 2 exp(Ha) increases 
exponentially with Ha 3> 1, where Ga is the Galilei number, and C c ~ 1 + sech(iJa) 
decreases from its non-MHD value of twice the free-su rface steady-state velocity to unity. 

As is also the case in channel flow ()Takashimalll996l ) , we recorded little variation in the 
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hard mode's critical parameters between small- Pm problems with insulating boundary 
conditions and the corresponding inductionlcss flows. On the other hand, we encountered 
considerable differences in the stability properties of the soft mode, manifested in the 
structure of eigenvalue contours in the (Re, a) plane, as well as the dependence of its 
critical parameters on Ha and Pm. Specifically, in problems with an insulating wall, our 
numerical results, supported by large-wavelength asymptotics, indicate that the a c = 
axis ceases to be part of the soft mode's neutral-stability curve, and the exponential 
growth of the critical Reynolds number becomes, for sufficiently large Ha, suppressed to 
a sublinearly increasing function. When conducting boundary conditions are imposed, 
Re c ~ (Ga/ PmY^Ha -1 becomes a decreasing function of the Hartmann number. 

The observed Pm-sensitivity of the soft instability was attributed to the strong-field 
behaviour of the participating inductionlcss mode (here called mode F), which, even 
though stabilised by the magnetic field, approaches neutral stability as Ha grows, and 
its energy content becomes almost exclusively gravitational. In particular, its decay rate 
and kinetic energy respectively decrease like Ha~ 2 and Ha~ A , where the latter scaling is 
consistent with a work balance between gravitational and Lorentz forces. The resulting 
near equilibrium is distinct from the quadratically increasing Lorentz damping experi- 
e nced by the s hear modes in the A, P, and S families (labelled according to the convention 
of Macklll976l ). In particular, it renders mode F susceptible to effects associated with mag- 
netic field perturbations, which are neglected in the inductionless limit, even when the 
magnetic diffusivity of the working fluid is large. 

Our analysis has identified two ways that nonzero magnetic field perturbations influ- 
ence the soft instability, both of which depend strongly on the wall boundary conditions. 
The first is through the component RePmJ X b of the Lorentz force associated with the 
steady-state current J and the perturbed magnetic field b. That force, which vanishes 
in the inductionless limit, results in a positive net energy transfer to mode F, leading 
in turn to the observed deviation of the critical Reynolds number from its behaviour 
in the inductionless limit. The fact that J depends, through the boundary conditions, 
on the wall conductivity, accounts for the different Re c results between insulating and 
conducting- wall problems. In particular, when conducting boundary conditions are en- 
forced, the magnitude of J increases without bound with Ha, and the resulting energy 
transfer to mode F causes Re c (Ha) to become a decreasing function. On the other hand, 
in insulating-wall problems J becomes constant throughout the inner part of the fluid 
domain, which is consistent with the comparatively milder modification of the soft mode's 
critical parameters. In this case, however, the boundary conditions are compatible with a 
stable, large- wavelength mode, which, as we have confirmed by means of large-wavelength 
approximations, is singular in the inductionless limit Pm \ 0. When Pm is nonzero this 
magnetic mode couples with mode F, causing the growth rate of the latter to become 
negative for sufficiently small a, irrespective of the value of the Reynolds number. The 
resulting large- wavelength instability suppression for all Re is the second major influence 
of nonzero magnetic field perturbations on the soft instability. As with the effects associ- 
ated with J, it too depends strongly on the boundary conditions. That is, the magnetic 
mode is absent from the spectrum when conducting boundary conditions are imposed, 
and in the same manner as inductionless flow, (for fixed Re > Ret) the soft instability 
takes place for arbitrarily small wavenumbers. 

Besides the large-wavelength magnetic mode, the spectrum of free-surface Hartmann 
flow with an insulating wall was found to contain a pair of travelling Alfven waves, char- 
acterised by a near-equipartition of the modal energy between the kinetic and magnetic 
degrees of freedom. At sufficiently high Alfven numbers, the upstream-propagating wave 
undergoes an instability where both Reynolds and Maxwell stresses are positive. Frozen-in 
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analo gues of the travelling Alfven modes (in the sense described bv lBetchov fc Criminate 



19671 ) were encountered in channel Hartmann flow with insulating walls, but are absent 



from conducting wall-problems in both free-surface and channel geometries. 

To conclude, the analysis presented in this paper highlights the important role played 
by the boundary conditions in free-surface Hartmann flow, and identifies potential short- 
comings of the inductionless approximation. Future work will explore signatures of the 
differences between the linear-stability properties of inductionless and low- Pm flows in 
fully nonlinear time-dependent simulations. 
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Appendix A. Large-wavelength approximations 

A.l. Non-MHD and inductionless problems 

A. 1.1. Formulation 

In non-MHD and inductionless free-surface problems, the Reynolds number Ret, at the 
bifurcation point (Reb,0) of the soft mode's neutral-stability curve, as well as the corre- 
sponding phase velocity C&, c an b e evaluated in cl osed form using regula r perturbation 
theory about a = 0. Following lYihl (|l963l . Il969h and lSmith fc Davisl (|l982h . we start with 
the expansions u(z) — u Q (z) + au\(z) + a 2 u 2 (z) + 0(a 3 ), a = a + aa x + a 2 a 2 + 0(a 3 ), 
and 7 = 70 + aryi + a 2 -f 2 + 0(a 3 ), which, when substituted into the modified OS equa- 
tion (|2~46|) . lead to a series of ordinary differential equations of the form 



D u r , 



D Ur. 



(Al) 



where n = 0, 1,2, . . . and := -(H 2 + Rej ) 
and depend on the solutions up to order n — 
orders we have 



Here the source terms s n vanish for n = 0, 
1 for n ^ 1. In particular, for the first two 



si := i?e(7i + iU) D 2 u - LRe(D 2 U)u + 2H X H Z D u , 
s 2 := i?e(7! + ill) D 2 u x - LRe(D 2 XJ)u x + <±H X H Z D u x H 
- (H 2 X + Re l0 ) Ul . 



{Re l2 + 2) D 2 u a (A 2) 



Similarly, the 5 —: Q boundary conditions (|2.48l a-c) and (|2.52p lead to 

«„(-l)=Du n (-l)=0, D 2 u„(0) = ^, (A3 a, 6) 

D 3 u„(0)^ 2 Du„(0) = Sl 4 l, u„(0) - 7oa« = Sg\ (A3 c,d) 
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where S, 



(3) 



vanish, and 



S[ 3) :=iD 2 J7(0)a , S. 



(3) 



iD 2 f/(0)ai -tto(0), 



S^ 4) := Ee(7i+if7(0))Du (0) +i{H x H z - ReDU(0))u (0), 



? (4) 



Re^ +iU(0))D Ul (0) +i(H x H : 
■f (3 + i?e 72 )DMo(0) 



.ReD^O))^) 



(A 4) 



GaRe ao, 



si 



(5) 



( 7 i+if/(0))a o , 



4 5) 



(71 + if/(0))ai + a 72- 



We express the general solution to (|A 1[) as 



i=0 



,(p) 



(A 5) 



-/i 2 DV 



0,^ 



where 4 are four linearly independent functions satisfying D 4 u\ - 
are constants, and Un(z) are particular solutions associated with the source terms s n . 
The constants K n i, as well as the expansion coefficients a n for the free-surface oscillation 
amplitude, are to be determined by systems of algebraic equations of the form 



A w, 



(A 6) 



that follow by substituting for u n in the boundary conditions (|A 3[) using (|A 5[) . Here 



,{h) 



An := 



Dm 



D 2 Un^ ' 



^ '(-1) 

.CO/ 



V 4» 



(-1) 



D 



.CO 



(-1) 



£» a 4 h) (o) 

(D 3 -^ 2 )^) 

4^(0) 



\ 




70 



(A 7) 



is a 5 x 5 matrix acting on five-element column vectors w n := (K u q, . 
the column vectors t„ , which vanish for n = 0, are given by 



S 



(3) 

n 

-M) 



D 2 ^ p) (0) 



^-(D 3 -/i 2 D)ui p) (0) 



V 5, 



(5) 



» 



(0) 



, ^„3,an) T - Also, 



(A 8) 



/ 



for n ^ 1 . Note that /4o is generally a nonlinear function of 70 , while the source vectors 
t n axe linear functions of "f n . 

Assuming that Aq has a go-dimensional right nullspace, denoted by ker(/4o) (as dis- 
cussed below, 70 will be chosen such that ker(/4o) is non-trivial), the solution to (|A 6p 
can be expressed as 



w n = RA n r , 



,(p) 



(A 9) 



where Ra is a Q x g matrix whose columns form a basis for ker(/4o), iT„ is a go- 
dimensional column vector of free parameters, and Wn is a particular solution associated 
with the source term t n . Therefore, introducing the notation v n := (u„, a„) T and := 
(u n P \o) T , as well as the matrix 



M 



,C0 



.(h) 



(A 10) 



50 



D. Giannakis, R. Rosner and P. F. Fischer 



the solution for the velocity field and the free-surface oscillation amplitude at n-th order 
in perturbation theory becomes 

v n = MR Ao n n + Mw%> + (A 11) 

In what follows, we choose the homogeneous solutions 

Uq := 1, := z, := (1 — cos(/iz))//i 2 , Ug := (jj,z — sin(/iz))//i 3 , (A 12) 

all of which are well-behaved in the limit fj, — > 0, giving, upon substitution into (|A 7p . 

/ 1 -1 (1-cos^)/^ 2 (sin^i - /^)/(U 3 \ 

1 — sin fi/fi (1 — cos fi)/fi 2 

A = 1 (A 13) 

[i 2 1 

V 1 -7o / 

At zeroth-order, the homogeneous problem AqWq = has a non-trivial solution only 
if A has a non-trivial nullspace, or, equivalently, 

det(/4 ) = cos(^)7 = 0. (A 14) 

The equation above has two distinct classes of roots, given by 70 = and 70 = —(Ft 2 + 
(2n + l) 2 7r 2 /4), where n = 0, 1, 2, . . .. Among these, only the zero solution can poten- 
tially be connected to a large-wavelength unstable mode, since the eigenvalues associ- 
ated with the first class of roots approach zero from below as a \ 0. Setting there- 
fore 70 = equips Aq with a one-dimensional nullspace spanned by the column vector 
£ := (0, 0, 0, 0, 1) T . Thus, the parameter vectors II n become scalars, playing the role of 
normalisation constants and, through (|A we obtain vq — (uo,ao) = #o(0, 1). We 
remark that channel problems do not admit asymptotically neutral solutions as a \ 0; 
in this case det(/4o) = 4(sin(/i) — /1 cos(/i)) / fj, 4 tends to 4/3 in the limit 70 — > 0. 

At higher orders in a one has to find solutions to inhomogeneous systems of equations 
of the form (|A 6j) . Here we will outline an inductive procedure which, given a zeroth- 
order solution, can be applied to obtain solutions at successively higher orders in a, 
and can also be generalised to treat the coupled differential equations (12.43|) governing 
nonzero- Pm problems (see appendix IA.2[) . 

First, assume that the eigenvalue coefficients 70, ... , 7 n -i and the corresponding eigen- 
vectors Vo, . . . ,v n -i have been evaluated to some order n — 1, where n ^ 1. More- 
over, assume that {vi}^^ are linear and homogeneous functions of go free parameters 
IIn-l,i; ■ • • i n n _x >qo , i.e. 

Vi = D iin _iJ7„_i (A 15) 

for 2 x </o matrices Do !Tt -i, . . . , D„_i.„_i and a go-dimensional column vector iT n _i := 
(77„_i.i, . . . , 27„_i go ) T . Under these conditions, the particular solution u„ , the boundary- 
condition source terms Sn > • • • > Sn and, by construction, the elements of t n at order n 
are also homogeneous linear functions of {n n ^ij}f?_ 1 . That is, we can write 

t n =T n n n - 1 , (A 16) 

where T n is a Q x go matrix. 

In general, a solution w n to (|A 6[) will only exist if t n lies in the range of /Ap, denote d 



by ran(/4o). According to the fundamental theorem of linear algebra (e.g. IStranal2005l) . 



ran(/4o) is the orthogonal complement, in the sense of the Euclidean inner product, of 
the left nullspace of Aq (i.e. the nullspace of Aq , denoted by ker(/4 r )). Thus, a solution 
to (|A 6j) exists if and only if L\ t n = 0, where La is a Q x go matrix whose columns form 
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a basis of ker(Ag )[J It therefore follows from (|A 16j) that a solvability condition for (|A 6j) 
is that the go x go matrix A n := T n has a non-trivial nuUspace, i.e. det(A fl ) = 0. Since 
det(A„) is a polynomial in j n of degree no greater than go (recall that the elements of t n 
depend linearly on j n ), the latter equation yields up to go distinct solutions for the n-th 
order expansion coefficient 7„. 

Denoting the dimension of ker(/4„) corresponding to a given solution for j n by q n ^ go 
(the procedure can be repeated for each of the roots of det(A n ) — 0), we now express the 
parameter vector iT„_i as 

II n -i=R An n n , (A 17) 

where Ra„ is a go x q n matrix whose columns are basis vectors for ker(A„), and JT„ := 
(IJ n ,i, ■ ■ ■ ,7J„ j9ll ) T is an updated vector of free parameters in the solution. Upon sub- 
stitution into (|A 15j) . (|A 17|) leads potentially to a decrease in the number of degrees of 
freedom in the eigenfunctions of order up to n — 1, as well as in wif^ (through its de- 
pendence on i„), from g to g„. Moreover, since the particular solution to (|A 1[) also 
depends linearly and homogeneously on iT„, it is possible to recast (|A lip as 

v n = MR Ao n' n + D n fl n , (A 18) 

where TI' n is a (provisional) g -dimensional column vector of free parameters and D n is 
a 2 x g„ matrix such that b n fl n — Mw^ + i4? . Although (IA 18|) may contain up to 
go + <7n arbitrary constants, the part of Il' n that is parallel to n n -i can be set to zero, 
since its only effect would be a renormalisation of the lower-order solutions. We therefore 
require that R]± U' n vanishes, or, equivalently, 

n n = E n ii n , (A 19) 

where E n is a go x (go — q n ) matrix whose columns form a basis of the left nuUspace of Ra„ 
and II n is a (go — g„)-dimensional column vector. If go happens to equal unity, one can 
set Il' n equal to zero. Inserting (|A 19[) into (|A 18[) . we then obtain v n — D n TI n + D n TI n , 
where D n :— MRA E n , or 

v n = D n , n II n , (A 20) 

~ ~ T ~ T 

where D n n := (D„, D n ) and 7T„ := (TI n , tl n ) T are a 2 x go matrix and a go-dimensional 
column vector, respectively. The lower-order solutions can also be written in terms of 7T„ 
using 

Vi = Di, n II n , D i<n := (0 2x ( qo - qn ), D i>n _i/?A„), . = 0,...,n-l, (A21) 

where 2 x( go -q„) denotes the zero matrix of size 2 x (g — g„). 

The matrices Dq^, . . . , D n ^ n fully specify the perturbative expansion to n-th order, 
which, as assumed above, is a linear homogeneous function of go parameters. The proce- 
dure can be repeated for successively higher orders in a, but, since the solutions for the 
expansion coefficients rapidly become complicated beyond first order, it is convenient to 
implement it in a symbolic-mathematics programming language so as to minimise effort 
and error probability. 

f According to the fundamental theorem, the row rank and the column rank of any matrix 
are equal. Thus, since Ao is square, the dimensions of its left and right nullspaces are equal, 
which implies in turn that La has size Q x go- 
lf The column rank of the go x q n matrix Ra„ is q n by construction (its columns are lin- 
early independent vectors). Since the row rank of a matrix is equal to its column rank, it 
follows that its row space (i.e. ran(/?J re )) is g n -dimensional, and its left nuUspace, ker(/?J ji ), is 
(go — g„)-dimensional. 
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A. 1.2. Results at first and second order 

We now apply the procedure described in §A.l.ll to evaluate the first and second-order 
corrections to the zeroth-order solution, (70,^0,00) = (0,0, 11$). Setting U(z) equal to 
the Hartmann velocity profile (equation (|2.26l a) with Co = Ci = 0), and the streamwise 
component of the applied magnetic field equal to zero (i.e. H x — and H z — Ha) the 
first-order expansion terms for the eigenvalue and the velocity eigenfunction are found 
to be 

«7! = l + scch(Ha), Ui{z) = — i77 sech(Ha) sinh 2 (Ha(l + z)/2)/ sinh 2 (Ha/2), (A22) 

while the coefficient for the free-surface amplitude, 01, vanishes. We remark that since 
71 is purely imaginary, the a — axis is part of the neutral-stability curve = Im(c) = 
Re(7i) + aRe(7 2 )-l-0(a 2 ) hi the (Re, a) plane. Moreover, the leading-order phase velocity 
Co := £71 is a monotonically decreasing function of the Hartmann number, with Co \ 
1 as Ha — > 00. As for the eigenfunction u\(z), it varies exponentially with the flow- 
normal coordinate, and like the Hartmann velocity profile, it possesses a boundary layer of 
thickness 0(1/ Ha) near the no-slip wall. In the vanishing magnetic field limit (Ha \ 0), 
Co is twice the steady-state velocity at the fre e surface, and u\ reduces to the quadratic 
function u\(z) = —in^l+z) 2 , as computed bv lYih ( 19631 19691 ) for free-surface Poiseuille 
flow. 

Proceeding now with the second-order approximation, 

Re coth(iJa/2) sech 3 (Ha)(2Ha(2 + cosh(2iJa)) - 3 sinh(2ila)) 



2 Ha 2 (cosh(Ha) - 1) 

8 Ga sinh 2 (Ha/2) (Ha - tanh(iJ'a)) 



Ha 3 Re(cosh(Ha) - 1) 

is the leading-order coefficient for 7 with nonzero real part, and governs the modal 
stability in the limit a \ 0. In particular, setting 72 equal to zero and solving for Re 
leads to (|4.1l a). quoted in the main text for the Reynolds number Reb at the bifurcation 
point (Reb, 0). That is, 72 is negative for < Re < Reb but becomes positive for all Re > 
Re b . For weak magnetic fields Re b /Ga 1/2 = (5/8) 1 / 2 + 191/(168 x lO 1 / 2 )^ 2 + 0(Ha 4 ) 
grows quadratically with Ha, but when the Hartmann number is large, Re b / Ga 1 ^ 2 ~ 
cxp(Ha)/ Ha 1 / 2 increases exponentially. Moreover, since 71 is independent of Re, the 
phase velocity C at the bifurcation point is equal to the zeroth-order phase velocity Co, 
in accordance with (|4.1l o) in the main text. We remark that the result Reb = (5Ga/ 8) 1//2 
for ze ro magnetic field strength is consistent with equation (38) in the paper by lYihl 
(1963), under the proviso that Re is replaced by Re' := 2Re/3 (Yih chooses the mean 



steady-state velocity as the characteristic velocity for reduction to non-dimensional form) , 
and the inclination angle 9 is substituted for Ga using (|2.4ip . 

In order to assess the relative importance of the formation of the Hartmann velocity 
profile versus the Lorentz force in the behaviour of Reb and C, we have carried out 
similar large-wavelength calculations for (physically unrealistic) problems with (i) the 
Hartmann velocity profile, but no Lorentz force (i.e. Ha set to zero in (|A 1I) - (|A4I) ). and 
(ii) the Lorentz force included, but the velocity profile set to the U(z) = 1 — z 2 Poiseuille 
solution. The results for the perturbation-expansion coefficients for the complex growth 
rate 7 are 



£71 = 1 + (Ha/2) 2 / smh 2 (Ha/2), 



(A 24a) 



72 
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Ga Re(432 - 24Ha 2 - 22Ha 4 + 5Ha 6 - 432 cosh(Ha) + 240Ha smh(Ha)) 

~ZRe~ 192ffa 2 sinh 4 (.ffa/2) ' 

(A 246) 



and 

«7i = 1 + 2(1 -sech(Ha))/Ha 2 , (A 25a) 

72 = -Ga(Ha-t&nh(Ha))/(ReHa 3 ) + Re(72 + 12 sech 2 (Ha) (4 + Ha tanh(Ha)) 
- sech( J ffa)(3(40 + 9Ha 2 ) + Ha(2Ha 2 - 3) t&nh(Ha)))/(6Ha 6 ), (A 256) 

respectively for cases (i) and (ii). As above, the phase velocity C& = Co = —171 at the 
bifurcation point follows directly from the first-order coefficients, while setting (|A 241 6) 
and (|A 251 6) to zero and solving for Re leads to the expressions 

R = 4ff a G a 1 / 2 (cosh(ffq) - 1) 

b ((18 + 5iJa 2 )(24-8 J ffa 2 + J ffa 4 )- 432cosh(iJa) + 240iJasinh(iJa)) 1 /2 1 > 

and 

Re b = 2Ha(3Gacosh(Ha)(Hacosh(Ha) - sinh(ffa))) 1/2 

x (-3(40 + 9Ha 2 + (40 + 9Ha 2 ) cosh(2iJa) - 12 cosh(3 J ffa) - 8Ha sinh(iJa)) 
+ cosh(iJa)(204 + 2Ha{3 - 2Ha 2 ) sinh(i7a)))" 1/2 , (A27) 

for the position of the bifurcation point on the a = axis. These two types of test 
problems have different strong-field behaviour, with 

Re b ~ Ga 1/2 (Ha/15) 1/2 exp{Ha/2), C b - 1 - (Ha/2) 2 exp(-iTa), 72 ~ -Ga/ZRe 

(A 28 a-c) 

and 

Re b ~ Ga 1/2 Ha 2 /3, C b -1~2/Ha 2 , l2 ~ -Ga/(ReHa 2 ), (A29a-c) 
respectively, for £/a 3> 1. 

A. 2. Nonzero-Pm problems 

The method described in ^A.l.ll can be used to study the large-wavelength limit of the 
coupled OS and induction equations (|2.43|) . with the addition that apart from the expan- 
sions for 7, u, and a, we write b(z) = bo(z) + ab%(z) + a 2 b2(z) + 0(a 3 ) for the magnetic 
field cigenfunction. Moreover, (|A 1[) are replaced by coupled differential equations, which, 
in the special case with flow-normal external magnetic field, have the form 

D 4 u» - Re la D 2 u n + HaPmT 1 ! 2 D 3 b n = s£\ 
D 2 b n - Rmjobn + HaPm 1/2 D u„ = s$ . 



(A 30) 



As in ^A.ll we are interested in perturbation order n ^ 2, where Sq^ = s^ — and 

s[ u) /Re := (71 + iU) D 2 u - i(D 2 U)u + i(D 2 B x )b - iB x D 2 6 , 
s^/Rm := (71 + iU)b - iB x u , 

sf /Re := (71 +iC/)D 2 wi - i(B 2 U) Ul + i(D 2 B x )b 1 -iB x D 2 b 1 (A 31) 

(72 + 2Re~ 1 ) D 2 u - 70 wo, 
s 2 b) /Rm := (71 + iU)b - LB a ui + (72 + i?m _1 )6 . 
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Among the boundary conditions, which are now 7 =: Q, the no-slip, shear-stress, and 
kinematic conditions, respectively (I2.48l a-c) have the same expansions as (|A 31 a. b, d), 
but the boundary condition for the normal stress (|2.48l c?) now yields 

D 3 w„(0) - Re l0 Du n {Q) + HaPm,- 1/2 Db„(0) = S£\ (A 32) 

where Sq 4 ^ = and 

S[ 4) /Re := ( 7l -i[/(0))D U o(0)-i(D[/(0))u (0)+i(DB !C (0))&o(0), 
S^/Re := ( 7l +iD'(0))D« 1 (0)-i(DD'(0))«i(0) + i(DB ie (0))(0)6i(0) (A 33) 
+ (3i?e _1 + 72) D u (0) + GaRe~ 2 a + HaPmRe^boiO). 

In problems with an insulating wall, the magnetic field boundary conditions (|2.50p lead 
in addition to 

D 6„(-l) = S< 6 > and D b n (0) = Sjp , (A 34 a, b) 

respectively, where sf ] = S [ Q 7) = and, for n e {1,2}, sf ] := & n _i(-l) and S^P := 
iD B x (Q)a n -i — 6„_i(0). When the wall is conducting, the boundary condition at z = — 1 
is replaced with b n (— 1) = 0, in accordance with (|2.5ip . 

Following an analogous approach as in (IA 5p , we express the general solution to (|A 30[) 
at order n as 

K,M =^iT ni (uf ) ,6f ) ) + (uW,6W), (A35) 

z=0 

where , bf 1 )}f = o are six linearly independent solutions to the homogeneous parts 

of (|A 30|) , and (vffl , bffi ) are particular solutions dependent on the source terms (s^ , Sn^ ) . 
Due to the high order of the differential equations involved, instead of seeking expres- 
sions for (nj ,&j for arbitrary 70, we shall set 70 = from the outset, which is the 
solution of interest for large- wavelength instabilities. We will verify subsequently that the 
resulting matrix Aq (the Pm > analogue of (|A 13|) ) possesses a non-trivial nullspace. 
Choosing then 

[u[ h \z)^\z)) := (1,0), (u[ h \z),b[ h \z)) := (z, -HaPm 1 ' 2 z 2 /2), 
(uf ) (z),bf\z)) := ((cosh(Haz) - l)/Ha 2 ,-Pm 1/2 (smh(Haz) - Haz)/Ha 2 ^ , 
(u^\z),bi h) (z)) := ( sinh ( ga z ) Z Ha z ; p m V2 cosh(gq g) - 1 - (Ha z) 2 /2 \ ^ 

{uf\z),bf r {z)) ~ (0,1), (ufW,^)) := (0, z) 

(A 36) 

as a set of linearly independent solutions to the homogeneous part of (|A 30[) . valid for 

70 = 0, and substituting (uo,6o) = Si=o -^"oiW ) m to the zcroth-order bound- 

ary conditions, leads to the homogeneous algebraic equations A w — 0, where w := 
(K 00 , . . . , Kq 5 , ao) T is a seven-element column vector, and A is a 7 x 7 matrix. In prob- 
lems with an insulating wall, Aq is given by 
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(A 37) 

where rows 1-7 respectively result from (|A3l a. b), (|A 32j) . (|A 341 a. 6), and (|A 3\ d). In 
conducting-wall problems, Aq has the same form as above but with the sixth row, origi- 
nating from the wall boundary condition for b, replaced by 



(Ay) 



HaPm 1/2 „ 1/2 sinh(ffa) - Ha n 1/2 2 + ffa 2 - 2 cosh(i?a) 




n Pm 1 - - Pm ' ■ - - 1 

2 ' Ha 2 ' 2i/ a 3 



The presence of all-zero columns in (|A 37[) . as well as in the matrix resulting from (|A 37|) 
with the substitution (|A 38|) . confirms our earlier assertion that setting 70 equal to zero 
endows Aq with a non-trivial nullspace. In the insulating-wall case that nullspace has 
dimension go := 2, and is spanned by the column vectors £ 1 := (0, 0, 0, 0, 0, 0, 1) T and 
€2 := (0) 0j 0, 0, 1, 0, 0) T , respectively corresponding to a free-surface displacement with 
zero velocity and magnetic field perturbations (as in inductionless problems), and a 
uniform flow-normal magnetic field perturbation with no change in either of u or a. In 
problems with a conducting wall the latter option is not available since the magnetic 
field must vanish at the wall; here ker(/4 ) is one-dimensional, and spanned by 

According to the discussion in §A.1.11 the two-dimensionality of ker(/4o) implies that 
in insulating-wall problems there exist up to two solutions for the first-order expansion 
coefficient 71, and carrying out the procedure outlined in that section establishes that 
this is indeed the case. However, unlike inductionless problems, the resulting expressions 
for 71 , which we denote by 7) and 7J , both possess negative real parts for all Ha > 0. 
Therefore, when Pm is nonzero all unstable inductionless modes acquire negative growth 
rate for sufficiently small a > 0, and the a = axis is no longer part of the neutral- 
stability curve Im(c) = 0. 

Explicit expressions for 7^ and 7^ ^ in terms of {Re, Ha, Pm} are complicated and 
not particularly illuminating. However, their salient properties are revealed by means of 
the series approximations 

(F ) 32Ha 2 Rm . / IZHa 2 l6Ha 2 \ ^ /rr4 , , A „„ , 

7f J = 0--1 2 T-)+0(Ha 4 ), (A 39 a) 

11 15(9 + 4i?m 2 ) V 90 5(9 + ARm 2 ) ) y h V ' 

( m) 2(3 + Ha 2 ) 32Ha 2 Rm , (2 TT 2 / 1 16 \\ ^ lTT 4 . 



3Rm 15(9 + ARm 2 ) \3 \90 5(9 + ARm 2 ) 



(A39&) 



and 



(F) _ Rm((cosh(AHa) - 1 - 16Ha 2 - 8Ha 2 cosh(2Ha)) tanh(ffa) + l&Ha?) 



li = 



6AHa 3 sinh 4 (iJa/2) cosh 2 (iJa) 
i(l + scch(Ha)) + 0{Pm 2 ), (A AO a) 
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(M) 2Hacoth(Ha) cosh(Ha) — 2Ha csch(Ha) + sech(Ha) ,. 

7i = -- 1 " " Utt \ V '- + 0(Pm), A40 6 

Rm cosh(Ha) — 1 

respectively valid for small Ha and Pm. Inspecting (|A 39[) and (|A 40|1 . we deduce that 
among the two solutions j[ M ^ is singular in the inductionless limit Pm \ 0, whereas 7 X 
tends to the result (|A 22[) for the first-order expansion coefficient for mode F obtained in 
appendix IA.1.21 For small Hartmann numbers the mode associated with 7 1 M '* , referred 
to in g3j] as the magnetic mode, has negative growth rate and propagates with phase 
velocity close to the (U) = 2/3 mean value of the Poiseuille profile. On the other hand, 

( F) 

mode F becomes neutral as Ha \ (that is, Re(7{ ; ) / 0), and propagates at twice the 
steady-state velocity at the free surface irrespective of the value of the magnetic Prandtl 
number. A similar procedure applied for problems with U = B = 0, but Ha ^ 0, yields 

7 { F) = 0, 7^ M) = -2H a coth(# a) /Rm, (A 41) 

indicating that the magnetic mode is stable even in the absence of a steady-state flow, 
whereas mode F becomes neutral to first order in a. 

(F) 

Because in insulating- wall problems Re^ ) is negative for all Pm > and Ha > 0, 
the second-order approximation in a is not relevant for the stability of mode F in the limit 
a \ 0, and for this reason we will not pursue it here. On the other hand, the analysis for 
problems with a conducting wall, where, as discussed above, ker(/4 ) is one-dimensional, 
leads to the result that, as in inductionless problems, the leading-order coefficient in the 
expansion for 7 with nonzero real part is 72, which, in addition to Re, Ga, and Ha, now 
depends on Pm. Solving for Re in the equation ^{Re, Ga, Ha, Pm) = 0, we obtain 

Re b = 8cosh(Ha) sinh(Ha/2) 2 Ga 1/2 (Hacosh(Ha) - sinh(Ha)) 1/2 

/(Ha(2(3 + (7 + 17Ha 2 )Pm) cosh(Ha) + ((13Ha 2 - 12) Pm - 6) cosh(3Fa) 

+ (Ha 2 -\)Pm cosh{5Ha) + Ha(2{& + IPm) sinh(Ha) + (4 + 5Pm) sinh(3i?a) 

- Pmsinh(5ifa)))) 1/2 (A 42) 

for the Reynolds number Reb at the bifurcation point of the neutral-stability curve, 
which in this case includes the a = axis. Moreover, the first-order coefficient 71 and, in 
turn, the phase velocity C& at the bifurcation point, are found to be given by the same 
expressions as in inductionless problems, namely (|A 22[) and (|4.1I &'1. respectively. For 
small Hartmann numbers, and provided that Pm is also small, Reb/ Ga 1 ^ 2 = (5/8) x l 2 + 
(191 - 25Pm)Ha 2 /(l68 x 10 1/2 ) + 0(Ha 4 ) increases quadratically with Ha, but for 
strong magnetic fields Re b ~ (Ga/ PmY^Ha" 1 becomes inversely proportional to the 
Hartmann number (cf. the exponentially increasing Reb in inductionless problems). 
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